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Si:CTION  1 
INTHOUrCTlON 


B.ACKGROl'Nl' 

A ma.ior  piV'bJem  durinp  the  desi^'n  of  a modern  hiph 
performance  tactical  missile  is  predict inp  and  modeling  the 
aerodynamic  forces  and  moments  it  will  experience.  This  prob- 
lem is  j)art  icula.r  ly  difficult  for  .t'.i.ssiles  wliich  must  execute 
extreme  manetivors  to  intercei't  hiph-g  targets.  The  mathematical 
model  of  the  airlrame's  aerodynamic  behavior  is  needed  for 
autopilot  and  guidance  law  design  and  for  a preliminary 
evaluation  of  the  missil<''s  perlormanci' . Thus,  the  initial 
design  effort  norntally  includes  many  hours  of  wind  tunnel 
testing  to  determine  an  aerodynamic  model.  Due  to  the' prob- 
lems of  tunnol  and  mount  interference,  dimensional  scaling, 
and  the  reduction  of  the  force  and  moment  data  to  a suitable 
analytical  form,  the  aerodynamic  model  obtained  may  be  highly 
inaccurate.  In  addition.  The  niissilo  airirame  configura- 
tion may  'oe  modi  lieu  during  the  design  to  improve  expected 
performance,  increase  control  effectiveness,  or  for  other  en- 
gineering reasons.  Modifications  to  the  aerodynamic  represen- 
tation accounting  for  such  changes  are  often  based  on  a limited 
amount  of  additional  wind  tunnel  testing  or  on  the  best  esti- 
mate's c'lf  the  af-rodynamic  designer. 

The  true  test  of  the  missile  simulat  loi'>  hkh  cl  i’St.d 
during  the  design  and  preliminary  evaluation  ci  r e.s  du'iiu, 
flight  test.  Toslflight  data  processing  involves  : vecon- 
struction  of  the  flight  using  data  taken  during  the  test, 
bis'ropanc  i es  between  exiiected  and  actual  llighi  behaviiu- 
in  many  cases  pc>int  to  inaccuracies  in  the  aerodynamic  model. 

In  addition,  unknown  effects  such  as  thrust  vector  misalignment, 


u f noini  na  1 thrust  levels,  and  llittht  ser.soi'  m-.ide]  iiiuert  a i n- 
tios  lead  t additional  d i screj’anc  ic'S . ^sin^:  the  ilifrht  test 
data  to  ir;iiiro\e  thi'  aerodynamic  iin'del  is  a crucial  jihase  of 
missile  desiKn  vhicii  t\i)jcally  results  in  airlramc  and/or 
autopilot  mod  i f 1 ca  t t ons  to  inijM'ovi.'  regions  o!  queest  iona  bK' 
per  forniunce . 

In  the  study  documen t ed  hero,  a method  ol  jiostfligiit 
data  ]>rocossing  is  invest  ig.at  ed  v.-hich  directly  addro.sscs  the 
airframe  modelinp  problem.  An  extended  Kalman  filter  is  do - 
Signed  to  estimate  a variety  of  aerodynamic  fiK)dt  1 and  moasuro^- 
ment  parameters  and  a structure  identification.  ]iroH’odiirc  is 
demonstrated  which  selects  among  candidato-  models.  A go^neral 
purpose  i.\-oiogrt-'fs-of-i  reedom  rnssiie  ilifht  test  model  is 
doivelojted  which  allows  flexibility  in  sons-.,'!’  selowticui  and 
parameterization,  as  well  as  rajnd  reprograntmi  ng  to  accomodate 
a variett  of  missile  airframes. 

The  concept  of  aerodynamic  model  determination  from 
flight  data  has  been  applied  with  success  to  aircr.ift  (Ref- 
erence 1).  An  aircraft  during  flight  ttist  can  be  well  in- 
strumented with  expensive  and  possibly  redundant  flight  data 
sensors,  which  can  be  used  repeatedly.  In  addition,  control 
surface  inputs  can  be  designed  w'hich  intentionally  separate 
complex  aerodynamic  ell.-cts  or  excite  aircraft  motion  in  a 
confined  flight  condition.  This  allows  relatively  simple 
aerodynamic  models  that  are  valid  for  special  cases.  Flight 
tests  can  bo  repealed  if  sensor  or  telemetry  problems  arise. 

By  comparison,  a missile  flight  test  program  pre- 
sents a more  difficult  aerodynamic  identification  problem. 

Because  ^^f  the  cost  and  the  destructive  nature  of  mis.sile 
flight  testing,  only  a relatively  few  flight  test.s  arc  scheduled. 
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Also  flight  tests  may  have  a variety  of  ob.lectives  in  addition 
to  verifying  cr  determining  the  aerodynamic  model.  Thus  the 
freedom  to  exercise  control  surfaces  for  the  purposes  of  aero- 
dynamic parameter  estimation  may  be  restricted.  Since  a single 
missile  flight  will  likely  involve  large  and  rapid  variation 
in  flight  condition,  aerodyna.mic  simplifications  applying  only 
to  restricted  flight  regions  must  be  minimized. 

The  ma:,or  difference  between  missile  and  aircraft 
aerodynamic  parameter  estimation  methods  lies  in  the  airframe 
modeling  detail  required.  Full  six-degrees-of ~f reedom  kinema- 
tics along  with  time-verying  thrust,  mass,  center  of  gravity 
and  inertias  are  modeled  in  order  to  characterize  the  tactical 
missile  throughout  the  operating  flight  regime.  This  modeling 
complexity  along  with  the  desired  ilexibility  and  generality 
of  the  resulting  postflight  data  processing  algorithm  places 
a considerable  burden  on  the  particular  parameter  estimation 
algorithm  used.  The  algorithm  must  be  computationi>.lly  effi- 
cient; it  must  produce  accurate  parameter  estimates  in  spite 
of  modeling  uncertainty,  and  user  convenience  must  be  maxi- 
mized. 


TECHNICAL  APPROACH 

The  application  of  modern  estimation  methods  to  the 
processing  of  measurements  from  complex  dynamic  systems  has 
become  commonplace  in  many  fields  of  engineering.  In  parti- 
cular, the  problem  of  system  identification,  or  parameter  esti- 
mation, has  received  much  attention.  Given  the  equations 
necessary  to  describe  the  dynamic  process  (the  state  equations), 
the  associated  measurement  equations,  and  the  system  and  mea- 
surement noise  statistics,  the  application  of  modern  estimation 
methods  is  well  understood  (Reference  2), 
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It  is  known  that  if  the  system  mathematical  descrip- 
tion is  accurate  and  if  the  equations  are  linear,  the  Kalman 
filter  gives  "optimal"  estimates  of  the  system  states.  For 
this  discussion, an  optimal  estimator  is  defined  as  one  pro- 
ducing lower  root  mean  square  trms)  estimation  errors  than 
any  oi her  estimation  technique.  In  addition  to  being  an  opti- 
mal estimator,  the  Kalman  filter  automatically  computes  the 
rrr,s  accuracy  of  its  estimates. 

The  linearity  restriction  mentioned  above  is  rarelv 
satisfied  in  realistic  physical  problems  such  as  the  one  con- 
sidered here.  Thus,  it  is  common  to  linearise  the  system 
equations  about  the  current  system  state  estimates  in  order 
to  dctcrm.ine  an  "appro.ximately  optimal"  system  estirnate.  If 
the  system  is  "nearly  linear"  in  the  region  which  contains 
the  system  states  and  estimates,  then  this  type  of  filter 
provides  a nearly  optimal  estimate.  This  approach,  used  here, 
has  been  highly  successful  in  a variety  of  engineering  appli- 
cations and  is  called  the  Extended  Kalman  Filter  (EKF). 

In  applying  the  EKF  to  parameter  identification,  the 
equations  of  motion  • — i.e.,  the  model  structure  --  are  assumed 
to  be  known  to  within  a set  of  uncertain  parameters.  The 
parameter  vector,  e.g.,  the  set  of  aerodynamic  coefficients,  is 
distinguished  from  the  vehicle's  dynamic  state  vector  --  posi- 
tions, velocities,  etc.  The  parameter  vector  is  adjoined  to 
the  dynamic  state  vector  to  form  an  augmented  state  vector. 

The  EKF  for  the  augmented  state  can  be  applied  in  a straight- 
forward manner.  After  measurement  processing,  that  portion  ol 
the  augmented  state  representing  the  parameter  vector  yields 
parameter  estimates.  Parameler  estimation  accuracy  is  obtained 
from  the  final  augmented  sta„e  covariance  matri.x. 
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A model  structure  identification  method  is  also  devel- 
oped and  evaluated  in  this  report.  The  method  requires  that 
several  candidate  structures  be  postulated  and  the  defining 
parameters  estimated.  The  system  models  are  statistically  com- 
pared using  information  stored  during  the  course  of  processing 
with  each  EKF,  and  the  most  probable  system  model  is  selected. 


ORGANIZATION  OF  THE  REPORT 

This  report  is  organized  as  follows:  Section  II  dis- 
cusses the  six-degrees-of-f reedom  airframe  dynamics  and  the 
aerodynamic  parameterization.  Frequently  encountered  missile 
flight  instrumentation  is  described  in  Section  ill  along  with 
the  associated  math  models.  Section  IV  describes  the  EKF,  pre- 
sents a partitioned  form  of  the  conventional  extended  Kalman 
iilter  tailored  specifically  for  the  large-scale  parameter 
identification  problem,  and  discusses  the  implementation  of  the 
EKF  for  the  airframe  identification  problem.  Section  V pre- 
sents aerodynamic  parameter  estimation  results  obtained  uti- 
lizing computcr-gcnorated  synthetic  measurements  from  a typi- 
cal short  range  air-to-air  missile.  A structure  identifica- 
tion method  is  developed  in  Section  Vi  and  test  cases  are  pre- 
sented which  demonstrate  its  ability  to  select  the  correct 
aerodynamic  model  structure  from  a family  of  candidate  models. 
Section  VII  discusses  an  application  of  the  identification  algo- 
rithm to  actual  flight  data  taken  from  an  aerodynamical ly  con- 
trolled airframe.  Section  VIII  summarizes  the  study  and  pre- 
sents the  major  conclusions. 


SECTION  II 

THE  MISSILE  AIRFRAME  MODEL 


OVERVIEW 


The  generalized  six-degrees-of-1 reedom  missile  air- 
frame model  used  in  this  effort  includes  12  fundamental  state 
variables : 

• Three  missile  inertial  velocity  com- 
ponents in  body  axes 

<s  Three  inertial  positions  to  locate 
the  missile  center  of  gravity  (C'G) 

e Three  rotation  rates  of  the  body  axes 
with  respect  to  inertial  space 

• Three  Euler  angles  to  locate  the  body 
axes  in  inertial  space. 

A generalized  aerodynamic  model  provides  force  and  moment  co- 
efficients as  a function  of  the  following  inputs' 

• Angle-of-attack 

• Angl e-of-sidesl ip 

• Mach  number 

• Control  deflections 

• Inertial  rotation  rates 

• Angle-of-a ttack  and  sideslip  rales 

• Total  relative  vel(.Mjity  maj^nitude. 

This  model  utilizes  a Taylor '.s  Scries  ex]  am  si  on  fcir  the  aero- 
dynamic lorcc  and  moment  coei  1 i ci ent -s  in  terms  ol  the  a.bove 


input  vur  i ubl  t.-s . 'J'unns  up  to  tliircl  ordt-r  an-  considered, 
with  approp.r  i :i  t e t orms  selected  liy  cunsiderinR  a variety  ol 
missile'  coni  ipurut  Ions  , us  well  us  utiliziin’;  the  syiruT'et  r i ca  i 
nature  of  th.r  gi  neral  missile  uirlruino  with  resjiect  to  a sin- 
gle jilano  of  sytaiiietry  (Reference  3).  The-  aorodynamu:  moments 
are  assuiiK'd  to  be-  reii'i'enced  to  a fixed  point  on  tho  missile 
uirlranu-,  and  must  be  translated  to  ucc'ouiU  for  CG  tiuvel 
wliile  thrusting.  The  llif.ht  i ime-h  i st  or  i os  of  three-  composite 
control  surlace  d'-iioctions  are  provided  --  the  pitch,  ya-^  , 
and  roll  contrei  tic  I 1 ec  t i ons  . Tiie  lt>vt'l  ol  mcideling  tiei  a 1 1 
assumed  is  typical  t)l  that  used  in  po^•;t  1 1 i glit  trajectory  rc- 
const  rut:  t i t'li  1 rtiit.  flir.ht  data. 

'riu-  o;irt.  h is  assumed  flat  and  m.-nrotating  ior  the 
duration  ul  t lit  missile  flight.  An  a tmt'Sj'ln-re  mtuiel  provide.s 
density,  pre.s'.,ure,  sju-od  oi  si'und  and  inertial  wind  velticitv 
a.s  a lunction  ol  altiiiuic.  Ftir  missile  i light  tests  at  Inf.li 
altitudes  the  wind  can  cont  r i luit  e si  gni  f icaiu  1 y tc'  tin-  (-llec- 
tivc  anr.  le-t)l --a  i t ack  . Ttic  t;lleeis  oJ  a t mos)iher  i c wind  on 
missile  aerody  nanii  c lorces  and  moments  is  accurately  i''])ri-- 
sent  ed . 

Th<-  missjlt  alrliamc  ma\  \n  ttii'ustinp  with  variabK- 
mass,  mniii-ni.s  tO  inertia,  and  CG  looalioio'.  Thrusi  valut's 
are  d<  t t-ri'.,i  m ti  l-y  an  iiijuit  table  oJ  sea  level  tliiust  Versus 
iiigln-time,  Willi  correci  ions  for  thru.io  varialion  w i i h .t  i t i - 
i udi  . 'i'lie  mass,  inertia  and  CG  proinrtie:  ai'c  cvii:iiun<\l 
j'orti'Uial  lo  tin  integral  ol  thrinst  t impulse),  ri  (g,!  .i 

injuits  v'lily  the  launch  and  hurneul  values  ■>'  m-  . :ii-  'i  . .i  . 
aiui  CG  local  i tin  . 


Tin  s 1 X -deg  rees-t>l  - 1 I'ot  dom  uiissile  imsl'  1 tKscrib  d 

^ieCl  !lld  el' til  r t 1 1 e.  l ■- 


aiiiiVf  is  lu't  exact  in  th.ii  a variety  el 
are  neglected  sinh  as  Innly  Viendiiig.  etc. 


1 1 j 


Cl  'liimi  lu  ) I r ' t ' e \ 


in  \ he  engir^eerinp  application  of  Kalman  filtering  methods  to 
account  for  these  known  (and  sometimes  unknown)  model  details 
by  assuming  random  inputs  (process  noise)  to  the  dynamic  equa- 
tions of  motion.  To  this  end,  six  additional  random  states  are 
added  t'  the  12  fundamental  airframe  states.  These  can  be  used 
to  rep’’i  sent  modeling  uncertainties  such  as  bending  effects, 
unsteady  aerodynamics,  wind  gusts,  or  control  deflection  mea- 
surement errors.  Three  random  force  state.s  and  three  random 
rrorr'ont  states  are  modeied.  The  rms  magnitude  and  bandwidth 
for  each  disturbance  state  are  required  inputs  to  describe 
the.se  randoia  effects. 


COORD IN.^TL  .hNP  NOTATION 

Two  coordinate  frames  are  used  frequently  in  the 
missile  airframe  model;  these  are  an  inertial  frame  and  a 
missile  body-fixed  frame.  The  direction  cosine  matrix  com- 
puted from  an  appropriate  set  of  Euler  angles  provides  the 
necessary  transformation  between  frames.  Superscripts  1 
(inertial  frame)  and  13  (body  frame)  are  used  to  indicate  the 
particular  frame  associated  with  the  coordinates  of  a given 
vector , 


The  inertial  coordinate  frame  assumes  a flat  non- 
rotating earth  and  is  fixed  at  launch.  The  origin  of  the 
l-fr;uTic  i:  ]oc;  led  on  the  surface  of  the  earth  below  the 
launch  point  (IJgure  1).  The  x and  y inertial  directions 
are  horizo  ital  whl.ee  the  z axis  is  down.  The  x axis  is 
oriented  in  the  dowiirange  direction  while  y is  crossrange 
to  the  right.  Noti  that  altitude  is  in  the  negative  z di- 
rection in  this  set  of  coordinates. 


EARTH’S  SURFACE 


Figure  1.  Inertial  and  Body  Coordinate  Systems 


The  body  coordinate  frame  is  located  along  principal 
body  axes;  i.e.,  all  cross  products  of  inertia  arc  zero  (Fig- 
ure 1).  The  x-axis  is  along  the  centerline  positive  forward. 
The  y-  and  z-axis  define  the  missile  control  planes;  i.e., 
missile  control  commands  are  given  in  the  x-z  plare  (j'itch 
plane)  and  the  x-y  plane  (yaw  plane).  The  origin  oJ  In 
body  frame  is  the  missile  CG. 


The  transformation  between  inertial  and  body  axes 
is  accomplished  via  the  direction  cosine  transformation. 

This  transformation  is  defined  in  terms  of  three  Euler  angles, 
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c = roll  unplo 
u = piicli  atiitudy 
; - yaw  uttitudo 

eci:r.;- 1 • • , descr  ipi  i on  ol  those  anp,los  is  contained  in  Rcl- 
orvirr  1.  The  resulting  direction  cosine  transformation  1 rom 
bod\  iM  inertial  ( ")  is  given  by  Table  1. 

'!  An  l:  1 . I'l  RECT  T ON  COS  I Nl.  TRANSFORMAT  I ON  MATE  I X ELEMENTS 


^ . j 

i t 

1 , 1 

cosc 

cost 

i,:i 

sin; 

s i n t' 

cos,  - cose 

si  n ; 

l , 3 

cos; 

sin? 

cosv'  + sin; 

Sint 

2 . 1 

COS0 

sin; 

2.2 

•s  in; 

sin? 

sine  ->  cos;- 

cos  e 

2,3 

cos-r 

sine 

sin^  - sine 

cos ; 

3,1 

- sinr* 

3,2 

sin?' 

cose 

3,3 

.1  ,ii„.  1 

cosc 

Ct-rtain  notation  conventions  used  in  the  missile  model 
r.r^.  ; a'esi'r. t e J uclow: 


• All  vectors  are  underscored 

• A dot  o\-cr  a quantity  indicates  its 
total  t i me -dcr i vat i ve 

• Elements  of  a vector  art?  subscripted 
by  the  appropriate  index,  i.e. . W3  is 
the  third  component  of  The  vector 


A "x"  between  two  three-element  vectors 
represents  the  cross  product,  i.e., 

w X yB. 

The  cross  product  equivalent  matrix 
w is  given  by 


I R - B 

-w-  I such  that  w x = wv 


The  quantity  g will  be  used  to  repre- 
sent the  gravity  vector  in  inertial  coor- 
dinates , i.e., 

I 

£ = 0 
g 


The  direction  cosine  transformation 
relating  the  I-  and  B-franies  is  given 
by  where,  for  example 

yB  = yl 


Thus  Cy  transforms  vectors  from  the  I-  to  the  B-frame.  The 

^ B I 

inverse  of  Cj  will  be  designated  Cg  so  that 

1 ~ Cfi  V 


MISSILE  EQUATIONS  OF  MOTION 

The  nonlinear  differential  equations  which,  must  be 
propagated  to  determine  the  12  f undP-tnental  missi.’c  :-?tates 
consist  of  the  four  interrelated  three- element  vector  equa- 
tions summarized  below: 


QACj  Ip 

m m 


v^ 
B - 


c?  s' 


w X 


w 


QAD  - OAc£  > Cj  + 


/ 


+ ^sini}  \v,.  + cos(}'  \v^)  tan0 

1 jCt  O 

cos(|'  - sin4> 

O 

\ W'o  + cosO  w,)/cos6  / 


+ 

— m 


(3) 

t4) 


Titj  sidof,  of  thes;e  lour  equations  represent  the  deriva- 

tives of  the  12  f undaiTiental  missile  states,  defined  as  follows 


« mi.ssile  velocity  with  respect  to  iner- 
tial space  in  body  axes 

r/  missile  CG  position  vector  in  inertial 
space 

w « missile  rotation  rates  with  respect  to 

~ inertial  space  always  understood  to  be 
expressed  in  body  coordinates 

V « a vector  made  up  of  the  three  Euler 

anples  which  locate  the  body  axes  orienta- 
tion with  respect  to  inertial  space. 


The  following  symbols  are  used  above  and  will  be  further  de- 
fined in  subsequent  subsections: 


0 = dynamic  pressure 

A ~ aerodynamic  reference  area 
Cj  aerodynamic  force  coefficients 
m - mass 

fij,  - thrust  vector  in  body  axes 

r^.  = acceleration  perturbation  input 

I.,  inertia  matrix  contuininc  inertias  I ,I  ,1 
N X V z 
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D = aerodynamic  reference  length 

C = uerodvnamic  moment  coefficients 
—m 

- \'ector  from  aerodynamic  reference 
point  to  the  CG 

m^^,  = thrust  moment  vector 

r = perturbation  moment  vector 

The  right  side  of  Equation  1 includes  the  summation  of 
all  forces  acting  on  the  missile.  The  five  teriiis  on  the  right 
side  of  the  equation,  from  left  to  right,  represent  the  follow- 
ing: 


• Aerodynamic  forces  in  body  axes 

• Thrust  forces  in  body  axes 

• Projection  of  gravity  onto  body  axes 

• Coriolis  contributions 

• Random  perturbation  effects. 

Equation  2 defines  the  rate  of  change  of  missile 
inertial  position  as  the  transformation  of  the  body-referenced 

r? 

velocity  v . The  right  side  of  Equation  3 represents  the  sum- 
mation of  all  moment.s  acting  on  the  airframe.  The  quantity 

is  the  diagonal  inertia  matrix;  the  five  terms  included 
within  the  braces  in  Equation  3,  from  left  to  right,  represent 
the  following: 

• Aerodynamic  moments  about  the  body-fixed 

aerodynamic  reference  point 

• Aerodynamic  moment  translali<-)n  to  CG 

1 ocat ion 

• Thrust  iiionient  about  the  CG 


Gyroscopic  coupling  terms 
Random  moment  perturbation  terms. 


Equation  4 represents  the  dynaniic.s  of  the  Euler  anple  states; 
Tlie.se  equations  are  derived  in  Reference  4. 


THRUST  AND  MASS  MODELS 

The  analytical  expressions  for  aerod^•namic  and  thrust 
forces  and  moments  are  in  some  cases  hiphly  airframe-specific. 
Thus,  the  general  purpose  airframe  model  proposed  here  may  not 
be  applicable  in  all  cases.  Nevertheless,  there  is  const der- 
able  comntonal i ty  between  the  majority  of  missile  airframe 
models  so  that  few  changes  will  be  necessary  to  represent  u 
variety  of  airframe  types. 

Variations  among  different  missile  airframes  are 
separated  into  two  categories, 

• Thrust -mass  characteristics 

• Aerodynamic  representation. 

The  time-varying  thrust  and  mass  characteristic.-.'  of 
an  arbitrary  missile  are  fundamentally  i elated.  The  nature  of 
this  relationship  is  used  here  to  minimize  the  inputs  required 
to  define  a specific  missile  configuration. 

A missile  thrust  versus  time  profile  at  sea  level  is 
assumed  in  tabular  form.  Linear  i nt  erj'ol  ul  i on  provides  a 
continuous  sea  level  thrust  time-history.  The  ma.s.s 


I i 


1 lov,  rail.'  is  iissuiiK'd  jirt.-poii  i ona  I so  that  muss  and 

mass-Ti.0  utL'd  i.  ha  na:;  l cri  st  i os  rati  be  nividoled  as  linear  luno- 
tiens  (.'!  the  thrust  intet^ral.  Tin.  missilv  tnas;s  eliarauter- 
ist  : CS'  requi  la  d are  ; 


Ma.ss  (m) 


Moments  ol  inert  ia  abc'Ut  body  axei 
(assuna'd  j-i' i lu- i i 'ul  axis)  1 , ] . I 

with,  diapv'nu.  inertia  tensor,  ■ 


Co  ii'dinates  ol  C(3  traw'l  Irom  the  aero- 
dynari'itc  J’Oler  n-'Ce  j'oint,  (c_S-a  three 
oicniont  voctoi  ). 


The  time-variation  o.l  tliese  stnen  n-.a ss~rel a t ed  vai’inhit's  is 
foaiputed  iron',  tlie  launch  (lull)  \alues,  the  burnriui  u'in]''ty') 
values  and  the  thrust  versus  tinn*  table.  In  particular, 


ml  t ) 


iHi  - 


LL^siilLlL 

Jo  " 


whore 


m( t ) = missile  mass  at  time  t 


= time  at  burnout 


- sea  level  thrust  at  time  t 


=■  launch  mass 


--  burnout  mass 


The  moments  of  inertia  and  CG  travel  are  computed  in  a similar 
manner. 


Sea  level  thrust  is  corrected  to  account  for  altitude 
variations  using  the  exhaust  gas  exit  The  resulting 

thrust  at  altitude  T is  given  by 


where 

= sea  level  reference  pressure 
^ (2116  Ib/sq  ft) 

pCalt)  = actual  pressure  altitude  from  the 
atmosphere  niodel* 

The  primary  thrust  direction  is  along  the  positive 
body  x-a\is.  Thrust  vector  rnisallgnnent s and/or  thrust  vector 
control  produce  lateral  forces  and  thrust  moments.  The  model 
presented  here  accounts  for  thrust  misalignment  but  does  not 
explicitly  consider  thrust  vector  controls  (TVC).  Note  that 
TVC  and  thrust  misalignments  are  modeled  in  a similar  manner, 
so  that  the  framework  for  a TVC  model  is  included. 

Figure  2 shov/s  a planar  representation  of  the 
thrust  moment  effects.  Note  the  body  axes  are  assumed  fixed 
at  the  missile  CG  which  may  be  located  arbitarily  with  respect 
to  the  geonictric  missile  centerline.  The  point  of  thrust 
attachment  to  the  missile  airframe,  as  well  as  the  aerodynamic 
reference  point , arc  assumed  to  be  fixed  along  the  missile 
geometric  centerline.  The  vector  eg  gives  the  CG  travel  from 
the  stationary  aerodynamic  reference  point.  The  c£  vector 

The  1962  ARDC  reference  atmosphere  is  assumed  in  this  study. 
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Figure  2.  Planar  Thrust  Misalignment  and 
CO  Travel  Diagram 

points  towards  the  CG  and  is  specified  in  body  axes.  The 
thrust  attachment  point  along  the  geometric  centerline  s 
referenced  from  the  aerodynamic  reference  point.  The  vari- 
able TAP  indicates  the  distance  from  the  aerodynamic  reference 
point  to  the  thrust  attachment  point;  the  positive  direction 
is  indicated  in  Figure  2. 


Small  angular  thrust  vector  misalignments  c and  r 

y ^ 

are  considered  as  positive  thrust  rotations  about  the  y and 
^ body  axes.  The  thrust  vector  in  body  axes  is  given  by, 


A vector  from  the  CG  position  to  the  thrust  attachniont  point  is 
given  by, 


-TAP 


CS) 
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Tl>ij  ttirii^^t  liioniont  about  the  CG  ir  K-’veu  by  : 


-TAl' 


b) 


AKHODY'v  \'ITr  MODKl.JNG 


S j .'i-d-'^vt'es-o  1- } reodf'in  aorodynanii  c modi;-]s  arc.  in  aio.^t 
casot;  !iL<zl!li,  cohicli*  ypec  i l‘i  c,  ut  i 1 i zi  in:  a \ar.iiM,\  of  one-, 
two-,  and  t hiww'-d  Lnioiisi  onal  tablos,  as  wdl  as  TaylC'r  series 
oxpan  .s  1 on  s , t I'i  ccuiopiet  r K exjmnsions  and  u‘,  tier  analytic  lunc- 
Tic'iis,  lo  cin.  ract  t-'i'ize  the  wind  tunnel  data.  The  comi.,cn  link 
between  i ho  \ jvious  aerodynamic  models  is  in  iln.'  i-ct  <'!  in- 


• I v U VJ  V U l. 


be  rv'latcd  :c 


orod\  :ian;i  c 


icrcos  and  tnc'nien  i >■ . These  i'a  ri  ables  a re  do  1' i nc'd  in  this  .si.'ction. 


The  relati\'o  velocity  vectory  defines  the  mo\-e- 
tnent  ol  the  mls.siie  throiiph  the  surrounding  air  mass.  The 
w.ind  velocity  ( ) .i.s  dofinod  as  the  movement  of  the  suio'ound- 
ing  air  mass  with  respect  to  the  inertial  coordinate  system. 
Consoquf’nt  1 y . t he  lelai  ice  velt.)eity  is  given  liy  , 


\ 

-r 


O 0) 


'J’he  volov.'itx'  assumed  to  hu?  afunctuui  only  ol  altitude. 

i\nJ  \ars'.i.o  .a.  tiiiiic  table.s  must  bo  sujipliod  to  tlie  model. 

The  rola-’\e  velc-city  is  expressed  in  t(.'rms  of  the 
relative  veloc.ity  magnitude  (v^.^)  and  two  angles  which  define 
the  position  of  the  relative  velocity  vector  in  body  axe.s. 
Two  .sets  of  relaii\e  velocity  vecto.r  orientation  angles  ai’o 
used  freciuently  in  missile  models.  Angl  o-of -a  t i ack  (a)  and 
angl  e-e  I - sidosl  i ]i  (f)  arc  used  exelusively  lOjr  aircraft  and 


for  winged  missiles  with  a dominant  single  plane  of  symmetry. 
For  these  vehicles  sideslip  is  controlled  to  be  small  and  this 
fact  often  simplifies  the  mathematical  representation  of 
the  yaw  plane  aerodynamic  model. 

For  missiles  with  two  planes  of  symmetry,  e.g.,  cru- 
ciform missiles,  aerodynamics  aie  frequently  modeled  as  func- 
tions of  aerodynamic  roll  and  total  angle-of-attack  . 

This  permi  t s the  use  of  the  known  vehicle  syrrjnetry  to  simplify 
the  aerodynamic  representation  (References  5 and  6). 

Angle-of-at tack  and  S'ideslip  are  sometimes  required 
To  be  Euler  angles  relating  the  body  uxe.s  to  a "wind  axes" 
system  --  i.c.,  a coordinate  system  with  one  axis  along  the 
relative  velocity  direction.  In  this  case,  a and  f.  are  defined 
by ; 


where 


Q = tan  (v  /v  . 

r,/  r ) 

\*  .i. 

(11) 

p-  = .s  i li  ^ / V ) 

' m 

(12) 

m 


2 ,1/2 
V ) 

^3 


(13) 


The  angl e-oi -attack  and  sideslip  can  be  used  to 
characterize  the  aerodynamics  of  missiles  wit)i  both  single 
and  dual  plane  oi  symmetry.  This  more  general  1 orniula t ion 
is  adopted  here.  Nevertheless,  it  is  recognizeo  ..hat  in 
special  cases  the  <)  ^ ) representation  may  load  to  reduced 

order  parameterizat  ions  and  thus  more  efficient  jiaramexer 
i dent i f icat ion . 
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An  alternate  definition  of  sideslip  is  frequently 
used  in  missile  modeling  where  the  ’’wind  axes"  are  rarely 
required.  This  definition  is  used  here  and  is  given  by, 

6 “ tan"^  (V  /v  ) (14) 

This  sideslip  definition  does  not  result  in  an  Euler  angle 
set.  The  similarity  between  the  a and  6 computations  of  Equa- 
tions 11  and  14  allows  some  simplification  of  the  aerodynamic 
model  for  cruciform  missiles. 

A generalized  missile  model  in  a and  B cannot  assume 
6 small.  In  many  cases  a and  B may  be  equally  large.  Never- 
theless the  assumption  of  at  least  one  plane  of  symmetry  -- 
the  xz  plane  — will  be  used  to  simplify  the  aerodynamic  model. 

The  model  given  here  assumes  all  aerodynamic  forces 
and  moments  to  be  linear  in  mach  number  about  some  reference 
v.alue,  H^.  This  us.sumption  limits  the  region  of  applicability 
of  the  model  to  flight  segments  where  Mach  variations  are  small. 
The  speed  of  sound  is  required  for  the  computation  of 
rnach  number  M.  The  speed  of  sound  is  computed  as  a function 
of  missile  altitude.  Mach  number  is  given  by: 


Both  aerodynamic  forces  and  moments  are  assumed  proportional 
to  dynamic  pressure  Q. 


Q 


(1C.) 


where  che  atmospheric  density  p is  computed  as  a function  of 
altitude  from  the  atmosphere  model.  In  addition,  the  aero- 


dynamic  force  in  Equation  1 is  proportional  to  some  specified 
reference  area  A while  the  aerodynamic  moments  (Equation  3) 
are  proportional  to  the  product  of  reference  area  and  reference 
distance,  AD. 

In  addition  to  a , 3,  and  M the  aerodynamic  coefficients 
are  functions  of  the  aerodynamic  control  surface  deflections. 
Missile  control  surface  deflections  are  commonly  governed  by 
three  independent  commands.  These  commands  are  designed  to 
produce  torques  about  the  three  body  axes  thus  achieving  con- 
trol in  pitcJi,  yaw,  and  roll.  The  pitch,  yaw,  and  roll  con- 
trol deflections  are  required  by  the  generalized  model  pre- 
sented here. 

Measurement  of  control  surface  deflection  requires 
instrumentation  of  actuator  response  — a measurement  fre- 
quently not  required  for  operational  missile  configurations. 

In  some  cases  actuator  response  measurements  may  not  be 
available  or  may  be  impractical.  A more  easily  obtained  mea- 
surement is  the  control  surface  command  which  results  from 
the  autopilot  subsystem.  These  commands  are  related  to  the 
actual  deflections  through  the  actuator  dynamics.  Actuator 
dynamics  are  typically  w-ide  bandwidth  with  respect  to 
the  dominant  airframe  frequency  response,  and  thus  may  be 
ignored  in  some  cases.  If  control  surface  command  rather 
than  actual  actuator  response  is  used  to  model  the  airframe 
input,  the  effects  of  actuator  dynamic  lag  should  be  investi- 
gated. This  can  bo  complicated  since  actuator  lags  may  in- 
volve calculation  of  aerodynamic  hinge  moments,  detailed  motor 
modeling  and/or  hydraulic  servo  models.  Tim  model  presented 
here  assumes  that  measurements  of  actual  control  deflections 
are  available. 


Six  independent  variables  required  ior  aerodynamic 
modeling  have  been  described  thus  far  --  angle-of-attack , 
sideslip,  much  number,  and  three  control  deflections.  These 
quantities  define  the  static  behavior  of  the  aerodynamics, 
i.e.,  a vehicle  fixed  with  respect  to  the  airflow  will  e.x- 
perience  forces  and  moments  a.s  functions  of  these  quantities. 
Most  wind  tunnel  studies  provide  only  static  empirical  force 
and  moment  relationships.  Additional  dynamic  forces  and 
moments  arise  as  a result  of  the  t ime- variant  ion  of  the  air- 
frmae  with  respect  to  the  airflow. 

The  dynamic  aerodynamic  effects  are  often  computed 
analytically  using  theoretical  aerodynam.lc  concepts.  Dynamic 
wind  tunnel  testing  methods  are  sonietimes  used  (Reference  7). 
most  common  dynamic  terms  considered  are  the  aerodynamic  damp- 
ing moments.  In  addition,  magnus  terms  are  sometimes  added 
which  couple  the  pitch  and  yaw  aerodynamics  when  roll  rates 
exist.  I'he  additional  variables  required  to  represent  the 
dynamic  terms  are  missile  rotation  rates  --  p q and  r represent 
ing  rotation  rates  about  body  x-,y-  and  z~axes,  re.spect i vel  y - 
as  Weil  as  a and  R.  The  standard  normalizing  factor  D/2v 
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is  used  to  normalize  the  dynamic  aerodynamic  terms. 

A generalized  mathematical  model  in  the  iorm  cl  Taylor 
series  expansions  relates  tlie  32  independent  variables  to  the 
(i  aerod.,'nami  c forces  and  moments.  This  model  is  deiined  in 
tcrn.iS  ol  time-invariant  parameters  in  order  to  make  use  of 
liarameter  ident  i 1 icat  ion  methods.  The  complete  aerodynamic 
nodel  is  given  by  Tables  2 through  5. 
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TABLE  5.  AERODYNAMIC  ROLL  MOMENT  --  C, 
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SECTION  III 

MISSILE  FLIGHT  TEST  MEASUREMENTS 


Data  taken  during  the  course  ot  missile  flight  tei.Mt- 
ing  for  the  purpose  of  post  flight  analysis  include  groui'd- 
based  missile  position  tracking  data  and  telemetry  records 
from  signals  generated  onboard  the  missile.  The  potential 
measurements  available  and  their  mathematical  models  art- 
discussed  in  this  section.  During  data  processing  ail  me;, sure 
ments  pre  assumed  discrete  with  a constant  data.  rate. 

ACCELERATION  MEASUREMENTS 

Three  body-fixed  accelerometers  are  assumed  to  pro- 
vide specific  force  measurements  along  the  body  x-,  y-,  and  z- 
axes.  These  sensors  are  typically  located  near  the  nose  of 
the  missile  and  thus  also  sense  the  effects  of  body  rotation^:. 
They  are  likely  to  be  "control  quality”  accelerometers  sub- 
ject to  large  biases,  scale  factor  uncertainty,  and  noist. 

The  accelerometer  measurement  is  modeled  by  , ■ 

a = [l+SFE  1 1=:^  4 d + w>;(wxr^,)  + v;xr^  ( 

— m L ^ — — —S'  — — sr  —a  —a 

where 

a = measured  acceleration 
— m 

f ® aerodynamic  force  vector 
—a 

==  thrust  force  vector 
w = missile  rotation  rate 

r = sensor  location  with  respect  to  the 

aerodynamic  reference  point,  in  body  axes 
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SFE  = 
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accelerometer  measurement  noise 
accelerometer  bias  error 
unit  identity  matrix 

diagonal  matrix  containing  accelerometer 
scale  factor  errors 


RATE  GYRO  MEASUREMENTS 


Three  body  fixed  rate  gyros  are  assumed.  Rate  gyros 
used  in  missile  control  are  subject  to  bias  errors,  scale 
factor  errors  and  acceleration  sensitive  errors.  The  rate 
gyro  model  used  is  given  by 


w = fl+SFE  1 w + [aSE]  + b„  + (18) 

~m  L gj  - -g  -g  -^g 

where 

= rate  gyro  measurement 

SFE  = diagonal  matrix  containing  rate  gyro 
® scale  factor  errors 

b = gyro  bias  drift  error 

O 

[ASE]  = 3^3  matrix  of  acceleration  sensitive 
drift  errors 

a = acceleration  in  body  axes  at  gyro 
location 


L - accelerometer  white  measurement  noise 
-^g 


The  gyro  and  accelerometer  instrument  packages  are  located  at 
approxiniR!.  ely  the  same  location,  so  that 


U -- 


- . --t  B I 

— + wx(wxrg)  + w^r^ 


(19) 


ANGLE-OF-ATTACK  AND  SIDESLIP  MEASUREMENTS 

Knowledge  of  the  relationship  between  the  angle-of- 
attack  and  sideslip  and  the  missile  dynamic  response  is  the 
key  to  determining  the  analytical  model  for  the  aerodynamic 
forces  and  moments.  In  some  cases,  angle-of-attack  and  side- 
slip can  be  determined  implicitly  using  the  accelerometer  and 
gyro  sensors.  Alternately,  the  angle-of-attack  and  sideslip 
may  be  measured  explicitly  using  the  appropriate  sensors.  Ac- 
celerometers and  gyros  measure  missile  characteristics  with 
respect  to  inertial  space,  whereas  angle-of-attack  and  side- 
slip are  functions  of  the  relative  velocity.  If  significant 
wind  is  present,  the  accelerometers  and  gyros  may  not  provide 
sufficient  angle-of-attack  and  sideslip  accuracy. 

For  missile  flight  tests,  the  a-B  sensor  consists  of 
a boom  projecting  ahead  of  the  missile  into  the  free  airstream, 
beyond  possible  interference  with  the  detached  shock  wave. 

Two  wind  vanes  are  attached  to  the  boom  and  allowed  rotational 
freedom  in  orthogonal  planes  (Reference  8).  The  vanes  remain 
aligned  with  the  relative  airflow;  thus  the  vane  angles  mea- 
sured with  respect  to  the  fixed  boom  provide  an  indication  of 
the  angle-of-attack  and  sideslip  at  the  missile  nose.  Mea- 
surement lags  result  because  of  vane  inertias  and  rotational 
friction , 


Various  attempts  have  been  made  to  calibrate  a-B  sen- 
sors by  inflight  calibration  methods  (References  9 and  10). 
These  studies  indicate  poor  agreement  between  manufacturer's 
calibration  results  and  inflight  verification,  as  well  as 
l.arge  discrepancies  in  the  predicted  sensor  lags.  Conse- 
quently, effective  use  of  the  a-B  sensor  may  require  identi- 
fication of  a variety  cf  sensor  parameters. 
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The  a and  B required  within  the  aerodynamic  model 
discussed  in  Section  II  are  derived  from  the  relative  velocity 
at  the  missile  CG.  The  ct  and  B as  measured  at  the  missile 
nose  will  include  the  effects  of  missile  rotations  with  re- 
spect to  the  airstream.  Consider  the  following  example: 

Missile  velocity  = 1200  ft/sec 
Pitch  rotation  rate  = 4 rad/sec 
a-6  sensor  to  CG  distance  « 3 ft 
Angle-of-attack  at  CG  = 0 


A discrepancy  in  angle-of-attack  at  the  CG  and  angle-of-attack 
at  the  nose  of  0.01  radian  (0.56  degree)  results.  Failure  to 
accurately  model  this  effect  may  account  for  the  in-flight 
calibration  difficulties  encountered  by  the  referenced  studies. 


In  addition,  a-6  booms  must  be  calibrated  to  account  for  bend- 

ing  under  load  (Reference  7),  and  they  may  lead 
tions  of  the  uninstrumented  aerodynamics. 

to  perturba- 

Ba.sed  upon  the  above  consideration,  the 
can  be  used  to  approximate  the  a-B  sensor  output 

following  model 

w r„ 
y B 

(20) 

w ro 

% = e * B 

(21) 

“cx  “l  “I.  “l  • “b 

(22) 

“s  “b 

(23) 

“m  “ a*SFE^)aj_  + 

(24) 

e,„  “ + be  'a 

(25) 

where 


^a’ 

a , 6 
m m 

SFE^.SFEg 


= angJ e-of -attack  and  sideslip  at 
nose  boom 

= distance  from  CG  to  nose  boom 

» dynamic  response  of  a-B  sensor 

= natural  frequency  of  a-6  sensor 

= damping  ratio  of  a-B  sensor 

= a-B  sensor  measurements 

= a-B  sensor  scale  factor  errors 

= a-6  sensor  bias  errors 

= wide  bandwidth  noise  on  a-6  sensor 

iVic  a 3 ur  cms  n t 


RANGE  ANI>  ATTITUDE  SENSOR  MEIASUREMENTS 

Several  sensors  can  provide  missile  range  and  atti- 
tude information.  An  inertial  reference  unit  (IRU),  the  on- 
board seeker,  and  remote  radar  and  theodolite  ranging  are 
discussed  here. 


An  inertial  reference  unit  consists  of  a package  of 
precision  accelerometers  and  gyros  which  provide  missile  posi- 
tion and  attitude  with  respect  to  inertial  space.  Because  of  cost , 
only  the  longer  range  tactical  missiles  have  an  operational 
requirement  for  an  IRU.  Nevertheless  an  IRU  may  be  provided 
as  part  of  the  flight  test  instrumentation.  Position  and 
attitude  measurements  from  an  IRU  are  subject  to  low  frequency 
drift  as  a result  of  initial  alignment  errors  and  instrument 
biases.  The  IRU  position  drifts  may  be  considerably  larger 
than  the  radar  position  uncertainty.  By  correlating  the  radar 
range  data  v/ith  the  IRU  data,  postflight  IRU  calibration 


can  be  performed,  resulting  in  estimates  of  instrument  biases, 
scale  factor  errors,  instrument  nonorthogonalities,  and  ini- 
tial alignments  errors  — along  with  an  optimal  estimate  of 
the  missile  trajectory  and  attitude  time-history.  This  pro- 
cedure can  be  performed  using  modern  estimation  methods  . 
without  a detailed  missile  airframe  simulation.  The  incTu- 
sion  of  the  postflight  IRU  calibration  within  the  airframe 
identification  algorithm  is  not  required:  however  the  capabil- 
ity of  using  the  results  from  such  an  effort  is  provided. 


Missile  test  ranges  may  utilize  a variety  of  radars 
spaced  along  the  missile  ground  track  during  a test  flight. 

In  addition  phototheodolite  sensors  may  be  available  to  pro- 
vide highly  a<-cuiaxe  triangulation  ranging.  Piecing  together 
the  various  radar  and/or  theodolite  measurements  is  a post- 
flight smoothing  problem  which  must  include  detailed  radar  and 
theodolite  error  models:  however,  again  an  airframe  model 
(Reference  11)  is  not  required. 


The  missile  inertial  position  estimates  from  combina- 
tions of  radar,  theodolite  and  IRU  measurements  will  be  con- 
sidered as  potentially  available  for  the  airframe  identifica- 
tion problem.  The  measurements  are  modeled  as  follows: 


r^  = + C 

— m — — r 


(2G) 


where 


- measured  missile  position  vector 
~ uncorrelated  position  measurement  error 


Modern  f i IterJ  ng/smoothing  methods  frequently  used 
for  postf light  trajectory  generation  provide  not  only  trajec- 
tory estimates  but  position  error  statistics  (i.e.,  the  .sta- 
tistics of  well.  Such  data  will  interface'  directly 

with  the  method  used  in  this  study. 


As  indicated  above,  missile  altitude  can  be  obtained 
from  the  IRU.  Attitude  information  may  also  bo  inferred  from 
an  onboard  seeker.  During  a missile  guidance  flight  test  an 
operationally  required  onboard  seeker  may  be  used  to  track  a 
target.  In  some  case.s  the  test  target  is  fixed  to  the  ground. 
During  such  a test  the  seeker  attitude  in  body  coordinates  to- 
gether with  the  target  location  provides  a measure  of  missile 
attitude  in  inertial  space.  The  attitude  inform.ation  and 
uncertainties  from  this  measurement  are  highly  seeker  specific. 
Errors  due  to  seeker  gimbal  dynamics,  radome  aberration 
effects  and  onboard  signal  ;•  ssing  roust  be  investigated. 
These  effects  may  be  modeled  using  process  and/or  measure- 
ment noise. 

Missile  attitude  information  from  an  onbo?.rd  IRU  or 
seeker  is  modeled  here  by  an  Euler  angle  measurement, 


where 

V = Euler  angle  measurement 
— m 

- uncorrelated  angle  measurement  noise 


(27) 


This  model  assumes  that  IRU  and  seeker  data  reduction  is  per- 
formed separately,  with  the  resulting  attitude  and  attitude 
uncertainties  stored  for  v.se  with  the  identification  algoriiam. 


SECTION  IV 


THE  EXTENDED  KALMAN  FILTER  FOR 
PARAMETER  IDENTIFICATION 


The  extended  Kalman  filter  is  an  extension  of  optimal 
(minimum-variance')  linear  filtering  theory  •co  problems  which 
involve  nonlinear  dynamics  and  measurements,  e.g.,  the  aerody- 
namic coefficient  estimation  task.  This  section  provides  the 
background  needed  for  understanding  the  EKF  algorithm  developed 
in  this  stud;\‘.  The  discussion  assumes  a basic  familiarity 
with  rai.dom  variables  and  state-space  notat  on;  additional 
details  can  be  found  in  Reference  2. 

KALMAN  FILTER  EQUATIONS 

To  apply  Kalman  iiltering  theory  to  any  estimation 
problem,  it  is  necessary  to  derive  a linear,  first-order, 
vector  dif lorential  equation  which  models  the  manner  in  which 
the  sy.stem  states  interact  and  propagate  as  tunotiens  of 
time.  For  linear  systems,  this  equation  has  the  general  form"*" 

x(t)  = F(.L)xCt)  + G(t)w(t)  + u(t)  (28) 

w'UTc  x(t)  is  an  column  vector  representing  the  system 

state,  E\t)  is  an  nxe  dynamics  matrix  which  defines  the 


*Some  overlap  in  symbol  usage  exists  between  this  s^ection  and 
Sections  II  and  111.  The  notiu  ion  in  thi.s  section  is  stan- 
dardized in  Reiorence  2,  to  which  the  reader  can  refer  for 
further  details. 


interaction  of  the  state  vector  components,  and  w(t)  is  a p>'l 
column  vector  of  white  gaussian  noise  inputs  such  that”^ 

Etw(t)]  = 0,  Cov[w]  - E[w(t)w(T)'^  ] = Q(t)6(t-T}  (29) 

The  matrix  G(t)  is  an  nxp  distribution  matrix  which  indicates 
how  each  component  of  w(t)  affects  each  component  of  the  sys- 
tem state  derivative,  and  u(t)  is  a n>^l  column  vector  of  known 
system  inputs.  Note  that  the  F,  G,  and  Q matrices  may  be 
time-varying.  For  a missile  airframe  model,  the  elements  of 
the  state  vector  x will  typically  be  missile  positions  and 
velocities;  the  elements  of  w will  be  random  inpiits  such  as 
airstream  turbul ence;  and  the  elements  of  u will  be  known 
inputs  such  as  missile  control  surface  deflections. 

At  discrete  instants  of  time  tj^  it  is  assumed  that 
measurements  of  linear  combinations  of  the  state  variables 
are  miade.  The  equation  describing  this  measurement  process 
has  the  general  form 


-k 


H.  X.  + V, 
K-K  — K 


(30) 


where  is  a vector  of  r measured  quantities  at  time  tj^, 
is  an  r>:n  observation  matrix  describing  the  linear  combina- 
tions of  state  variables  which  comprise  in  the  absence  of 
noise,  and  v^  is  an  r-vecLor  of  zero  mean  gaussian  measure- 
ment errors  with  a covariance  matrix  R.  defined  by 


(31) 


♦The  symbol  E[  ] denotes  mathematical  exj^cctat  ion ; Covlw] 
denotes  the  covariance  matrix  of  w-. 


At  any  time  t,  the  objective  of  optimal  estimation 
theory  is  to  process  all  the  measurements  taken  up  to  that 
time  and  produce  an  estimate  x(t)  of  the  system  state  x(t) 
having  minimum  error,  in  a statistical  sense.  The  optimiza- 
tion ciiiorion  most  often  chosen  is  that  of  minimizing  the 
mean  square  estimation  error.  The  minimum  mean-square  error 
estimate  is  c-ilculated  with  the  Kalman  filtering  algorithm. 


As  measurements  become  available,  there  is  an  instan- 
taneous change  in  the  knowledge  of  the  state  Denoting 

the  optimum  estimate  of  x(tj^)  just  prior  to  the  availability 
of  as  '■'^nd  the  optimum  estimate  of  the  state  vector 

immediately  .ifter  processing  as  Xk''"*'^’  Kalman  filter 

generates  the  optimum  estimate  of  the  system  state  according 
to  the  following  .it  gorlthm : * 


x(t)  = F(t)x(t)  + u(t)  ; x(t^_^)  = 

ikV)  - iK<->  - 

where  Eq . 4.1-6  is  used  to  calculate  the  estimate  between 
measurements  and  Eq . 4.1-6  is  used  to  update  the  estimate 
when  new  data  is  received  at  each  time  t,  . 


(32) 

(33) 


The  quantity  In  brackets  in  Equation  33  is  Known  ns 
the  measurement  residual  r . This  vector  quantity  is  the 
diifeiencc  betweer’  the  actual  measurement  and  the  optimal 
estimate  of  the  n.eisureroent  ^k?-k^“^'  important  attribute 

of  the  Kalman  filter  is  the  whitening  of  the  residual  process. 
It  can  be  shown  that  a Kalma-n  filter  produces  a residual  proc- 
ess  such  that 


♦Only  the  continuous  form  of  the  Kalman  filter  with  discrete 
measurements  is  considered  here. 


0 


il]  ' 


In  addition  it  Is  known  that 


^2 

Cov  [r^  rj]  . [h^Pj,(-)hT  R^] 


The  "whiteness"  of  the  residual  process  and  the  consistency 
of  the  residual  magnitudes  with  the  predicted  residual  stand- 
ard derivation  is  an  excellent  test  for  proper  filter  behavior. 


The  nxr  matrix  K-  is  the  Kalman  gain  mat.'ix.  Let 

K 

x( t ) denote  the  error  made  in  estimatiug 

x(t ) = ^(t ) - x(t ) ( 34  ) 


and  let 


P(t) 

= Covj^xC 

t)]  * 

E^x(t  )x(t)'^J 

(35) 

The  gain  matrix 

is  t 

hen  computed  using  the  following 

equal: 

P(t) 

= F(t)P( 

t)  + 

P(t)F(t)'^  + G(t)Q(t)G(t)'^ 

(36) 

with  P(ij^_ 

and 

II 

(37) 

Pk(-) 

Vk] 

(38) 

where  J’(t)  just  before  the  measurement  at  time 

and  is  PCt)  just  after  Equation  36  is  used  to 

calculate  the  estimation  error  covariance  between  the  measure- 
ments; Equations  37  and  38  are  used  to  calculate  the  Kalman 


pain  matrix  for  use  in  Equation  33  and  to  update  the  estima- 
tion error  covariance  matrix  when  a measurement  arrives. 

Figure  3 illustrates  the  procedure  for  prc)cessing 
measurements  with  the  Kalman  filter.  The  algorithm  has  two 
distinct  phases.  Equations  32  and  36  describe  the  time  evolu- 
tion of  the  state  estimate  and  its  error  .statistics  between 
measurements  under  the  influence  of  system  dynamics  and  noise. 
This  process  is  commonly  referred  to  as  extrapolation.  Equa- 
tions 33.  37,  and  38  indicate  how  the  estimate  and  its  error 
covariance  are  updated  at  the  measurement  time  to  reflect  the 
new  information  available.  The  algorithm  is  optimum  in  the 
minimum  mean-square  error  sense  as  long  as  the  assumed  math- 
ematical model  for  the  system  is  accurate. 
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Figure  3. 


Measurement  Processing  with 
the  Kalman  Filter 


A unique  feature  of  the  Kalman  filter  is  that  a mea- 
sure of  the  estimation  accuracy  is  automatically  provided  by 
the  algorithm.  The  matrix  P(t)  is  a complete  description 
of  the  error  statistics  of  the  estimate.  In  paiticular,  the 
square  roots  of  the  diagonal  terms  of  P(t)  represent  the 


: 4 


standard  deviations 
components  of  x(t.).-. 
by  Equations  36,  37 
Kalman  gains  can  be 
ment  sequence  . 


of  the  error  obtained  in  estimating  the 
Note  that  P(t)  is  specified  for  all  time 
and  38.  Thus  estimation  accuracy  and  the 
precomputed  independently  of  the  mease  re- 


in sumrr.ary , the  following  conditions  must  be  met  to 
im])iement  an  optimum  Kalman  filter; 


• The  system  must  be  linear 

• The  matrices  F(t),  0(t),  and  u(t)  must 
be  known  functions  of  time 

• The  vector  input  w(t)  must  be  a gaussian 
white  noise  process  with  known  covariance 
matrix,  Q(t)6(t-T),  and  zero  mean 

• The  measurements  must  obey  Equation  30, 

and  H,  must  be  known  for  all  k 
k 

• The  measurement  errors  V],;  must  be  a 
gaussian  white  seq.uencs~with  covari- 
ance niatrix  and  zero  mean 

• To  initialize  the  filter  equations, 
the  Initial  statistics  of  x must  be 
known . 


If  the  missile  airframe  identification  problem  could 
be  put  into  a form  which  met  all  of  the  conditions  listed 
above,  the  design  of  an  optimal  estimator  would  involve  only 
the  direct  implementation  of  the  Kalman  filter  equations. 
However,  the  airframe  dynamic^  and  measurements  considered 
here  are  nonlinear;  the  linearity  requirement  is  violated. 
Furthermore,  the  objective  of  estimating  the  uncertain  aero- 
dynamic coefficients,  which  is  tantamount  to  e.stir.iating  pa- 
rameters of  the  matrix  F in  Equation  28,  introduces  additional 
nonlinearity.  One  means  ol  overcoming  these  problems  is  the 
extended  Kalman  filter  described  in  the  next  subsection. 


EXTENDED  KALMAN  FILTER  EQUATIONS 

Since  the  problem  under  consideration  cannot  be 
realistically  modeled  as  a linear  system,  a nonlinear  estima- 
tion technique  must  be  used.  The  method  selected  here  is  the 
extended  Kalman  filter.  The  latter  is  e.ssentially  a conven- 
tional Kalmar  filter  design  applied  to  a mathematical  model 
obtained  by  linearirting  the  system  about  the  current  state 
estimate.  The  structure  of'  this  algorithm  is  illustrated  in 
Figure  4.  Note  that,  because  of  the  linearization  proce- 
dure, the  covariance  calculation  is  now  dependent  upon  the 
state  estimate.  Consequently , it  is  not  possible  to  precom- 
pute the  cuva  iance  matrix  and  Kalman  gains,  as  functions  of 
time,  since  they  are  dependent  upon  the  measurement  data.  The 
extended  Kalman  filter  yields  very  nearly  optimal  estimates 
if  the  linearisation  is  accurate,  i.e.,  a.s  long  as  the  state 
estimate  is  close  to  the  true  system  state. 
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FiK'ure  4 


Structure  of  an  Extended  Kalman  Killer 


i i 


A reasonably  general  mathematical  model  lor  nonlinear 
stochastic  systems  is  given  by  the  equations 


x(t)  * f(x(t),  t]  + G(t)w(t) 


Sk  = ik[2"k>]  - ; k - 1,2. 


(39) 


(40) 


where  f and  hj^  arc  nonlinear  differentiable  functions  of  the 
state  vector  x.  and  w(t)  and  are  zero  mean,  independent 
gaussian  white  noise  processes  having  spectral  density  and 
covariance  matrices  Q(t)  and  . respectively.  The  measure- 
ments z.  are  taken  at  discrete  times  t,  . 


The  first  approach  one  might  use  to  derive  a filter- 
ing algorithm  for  x(t)  in  Equation  31)  is  to  linearize  the  non- 
linear functions  ^ and  h^  about  an  appropriate  known  refer- 
ence trajectory  x(t),  and  then  apply  conventional  linear 
estimation  theory,  i.e.,  the  Kalman  filter  discussed  in  the 
last  subsection.  Thus,  denoting  expressions 


f(x.  t)  £ f(x,  t)  + -j- 


(x-x) 


(41) 


X,  ®x 
— K — 


(42) 


may  lu;  substituted  into  Equations  39  and  4(>  in  d.n)  ' Ui 
corresponding  Kalman  filler  whi  eh  e.stimat  es  in.  in  x 

from  the  i-eicrence  trajectory.  When  the  rcj-.  i i trajecluti 
is  chu.sen  to  be  the  c.urrcni  best  estimate  ol  i t,i.  state  x(t), 
the  rosuliing  algorithm  is  known  us  an  oxt<,uucu  Kalman  ti1t('r 
(EKF);  the  mechanization  equations  for  the  latter  are  giia-n 
in  Tabic  G (.see  also  Reference  2). 


sc-yaMiv  oy  thk  extejided  kalman  filter 


• TABLE  6- 


Tbf  rriutrl/.  >’ 
proximal- ion  to  1. ho  o.s 
Xj. (-)  and  dv.uoi. 

lions  :il  l 11U‘  t,  J-'.'-L 
Note  111:  III  o,; n 

filter  outlined  in  th 


(i)  ill  Tahi'*  i>  as  a f i r£.t-oi'd'0'  ap- 
t.  in. at  ton  error  covari  a.n-i  e z'jatrjx,  and 
1'  the  .v.ol'Jtiors  to  the  prono-joit  ion  e^ua- 
In'L.ire  the  ineaaui'oM.ciit  is  ;o’ocessod. 

ati'.'iit  vediu.o  to  the  rpiMma,]  Jial'iiiui 
e Iasi  re,'cti'.Pi  if 


fixu).  t)  = F(,l)x(t)  -1  t ) (4d) 

and 


"k2'T> 


(4-1) 


1 MPLf:MENTATI ON  CONS  IDERAT I ONS 


THE  PARTITIONED  FORM 

The  application  of  the  EKF  to  parameter  identifica- 
tion is  straiKhtf orward . An  augmented  state  is  defined  as 


where  is  the  dynamic  state  vector  and  a is  the  uncertain 
parameter  vector.  The  required  first  order  differential 
equation  for  x is  given  by 


V s 

pf(Xj.(t),  a,  t)*! 

4. 

"G  w(t)' 

A 

[a . 

1 

io 

0 

For  large  state  vectors  the  primary  computational 
burden  of  the  EKF  is  in  the  propagation  of  the  state  covari- 
ance between  measurements  (Equation  36).  As  an  example,  con- 
sider an  airframe  parameter  estimation  problem  containing  50 
parameters  and  18  dynamic  states.  For  this  sample  case  the 

augmented  state  is  of  length  68.  The  augmented  state  covari- 

* 

ance  matrix  will  have  2346  unique  elements.  Thus  an  extended 
Kalman  filter  lor  a state  of  order  68  would  have  to  numerically 
integrate  2346  + 68  = 2414  scalar  equations. 

The  numerical  integration  burden  can  b>-  rc-duc  ('d  b^ 
partitioning  the  augmented  state  covariance  mairjx  into  a 


*An  nxn  symmetric  covariance  matrix  contains  n(n+l)/2  unique 
elements . 
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dynamic  state  covariance  . a parsuneter  covariance  P 


ss 


aa 


and  the  state-parameter  "cross  covariance"  Pg^-  The  aug- 
mented state  covariance  P is  given  by 


P = 


p 

' p 1 

1 P 

^ ss 

1 sa 

r p 

sa 

1 aa 

1 J 

The  resulting  equations  become 


! P + P + f P"^  + P„„  fT  + GQg'^ 

— >.g  ss  ss  — Ag  — £ sa  sa  —a  ^ 


(4C) 


P = f r + f P 
su  —X  c,'.  —a  aa 

-s  - 


(•■17) 


'“aa  ' >»' 


where  element5  oi  iiu'  matrices  f and  £ are  given  by, 


3fi(2Ls>  t) 


fial  . . 

t -Jlj 


?a" 

J 


Note  that  the  parameter  covariance  P is  constant  between 

aa 

measurements  and  thus  requires  no  integration.  For  the  50- 
parameter,  18-state  example  the  numerical  integration  require- 
ments are  now  1089  scalar  equations,  or  a savings  of  over  50 
percent  . 


I 


I 
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PARTIAL  DERIVATIVES 


Four  partial  derivative  matrices  are  required  by 
xhe  partitioned  foim  the  EKF. 

• Partials  of  dynamic  state  derivatives 
with  respect  to  dynamic  states  fjj 

« Partials  vof  dynamic  state  derivatives 

v'lth  respect  to  parameters 

• Partials  of  measurements  with  respect 

to  state  h„ 

— xr 

« Partials  of  measurements  with  respect 

to  parameters  h 

These  partial  derivatives  are  computed  analytically 
dvxring  the  course  of  the  numerical  integration.  Conventional 
though  tedious  algebraic  methods  are  used  in  the  derivation 
of  these  equations.  Because  of  tho  required  generality  of  the 
resulting  filtering  algorithm, no  simplifying  assumptions  were 
made . 


In  addition  to  the  analytic  partial  derivatives,  a 
numerical  partial  derivative  calculation  has  been  utilized. 
The  latter  is  useful  for  two  reasons.  First,  it  offers  an 
independent  verification  of  the  analytic  partial  derivatives, 
and  second  it  can  provide  an  indication  of  the  linearity  of 
the  dynamic  equations  in  a region  about  the  augmented  state 
estimate . 


The  numerical  partial  derivative  is  based  upon  a two 
sided  perturbation  of  the  function  to  be  differentiated  — i.e., 
both  plus  and  minus  perturbations  are  used.  The  perturbation 
step  size  is  determined  based  on  the  standard  deviation  of 
the  dynamic  state  or  parameter  with  respect  to  which  the  deri- 

I 

vative  is  to  be  taken.  Thus  linearitv  of  the  dynamics  in  a % 


statistically  defined  region  about  the  state  and  parameter 
estiiTiate  is  assumed  in  the  calculation  of  the  numerical  par- 
tial derivatives.  Recall  that  accurate  linearization  is 
one  of  the  fundamental  assumptions  required  in  order  for  the 
EKF  to  be  a near-optimal  estimator.  Thus  correlation  between 
the  analytical  and  numerical  partial  derivative  calculations 
provides  both  an  indication  of  the  correctness  of  the  analy- 
tical partials  as  well  as  providing  a necessary  condition  for 
optimality  of  the  EKF  parameter  estimates. 

IDENTIFICATION  OPTIONS 

In  the  most  general  from  of  the  missile  niOdel  pre- 
sented in  Sections  II  and  III  there  are-*  22  dynamic  states,  1*1 
measurement. s , and  approximately  150  parameters  to  be  estimated. 
A summary  of  the  quantities  contained  in  each  of  these  cate- 
gories is  given  in  Table  7.  The  filtering  algorithm  is 
structured  so  as  to  allow  subsets  of  the  most  general  case  to 
be  exercised  with  particular  attention  given  to  allowing  a 
flexible  parameter  set.  The  aerodynamic  parameters  are  im- 
plemented so  that  particular  parameters  may  be  selected  for 
estimation  while  the  remainder  may  be  fixed  at  preset  values. 

It  is  anticipated  that  less  than  50  parameters  will  be  esti- 
mated at  one  time  for  any  particular  case. 

Certain  options  are  available  which  reduce  the  com- 
putational burden  during  the  calculation  of  the  analytical 
partial  derivatives.  These  options  are  listed  in  Table  8. 

COMPUTAT  rONAI.  F.EQU  I REMENTS 

The  present  algorithm  is  coded  in  FORTRAM  IV  and 
implemented  in  single  precision  in  an  IBM  360/14,5  computer. 

Core  requirements  are  approximately  340K  IBM  360  bytes  based 





TABLi:  7,  DYNAMIC  STATES,  MEASUREMENTS  AND  PARAMETERS 


— 

1 

i 

1 

, 

; 

DESCRIPTION 

eynami.;:  St  it  os 

1-3 

Missile  Velocity  in  Body  Axes 

4-3 

Missile  Position  in  Inertial  Axes 

7-  C 

Massile  Body  Rotation  Rate 

10-12 

Missile  Euler  Angles 

13-15 

Missile  Acceleration  Perturbar  .on  Components 

1 3- 1 3 

Mis.sile  Moment  Peituroaiion  Components 

19-20 

a Sensor  Dynamic  States 

21-22  * 

6 Sensor  Dynamic  States 

Measureraems 

1-3 

Bodv  Fi.xed  Accelerometers 

4 - C 

Body  Fi.ved  Kate  Gyros 

7-8 

c.-  5 Sensor 

9-11 

Position  Measurements 

12-14 

Attitude  Measurements 

Potential  Paran*?tor 

Li  St 

1-4  0 

Aerodynamic  Force  Coefficients 

41-100 

•Aerodynamic  Moment  Coefficients 

lGl-103 

'.Vind  Bias  Uncertainty 

104 

Atrt.c-spner  ic  Density  Uncertaintv 

105-106  * 

Thrust  Misalignmeui 

107  * 

Thrust  Bias 

108  * 

Mass  Uncertainty 

109-111  ■'(> 

Moments,  of  Inertia 

1 12-114  * 

CG  Docarions 

115-117  * 

Perturbation  Force  Time  Constant 

118--120  * 

Perturbation  Moment  Time  Constant 

121-123  )0 

Pei' l urbat  Ion  Force  Scale  Factor 

124-126  * 

Perturbation  Moment  Scale  Factor 

127-129 

Acceleronim  er  Bias 

130-132 

Accelerometer  Scale  Factor 

1 33-135 

Gyro  Bias 

I"''  138 

Gyro  Scale  Factor 

1 3. .-146 

Gyro  G Sensitivities 

14  9-;.  50  * 

a-6  Sensor  Natural  Frequencies 

151-152  * 

a-5  Sensor  Damping  Ratios 

1 53-154 

a -8  Sensor  Bias 

15  5-15G 

S‘:,'nsor  Scale  Factor  Errors 

L.  . ...  ...  ..  ... 

•+T*.\ost-  opiious  !?rt’  lun  ]n't-.-«cnl  ly  imjilomcnted  in  the  idem  i 1 icul  ion 
n i.'.'ier  1 I hni , 


TT  iiwm — r^rr*  -rrnjiii  'i 


TABLE  8.  ALGORITHM  COMPUTATIONAL  OPTIONS 


OPTION 

EXPLANATION 

No  Wind 

All  wind  calculations  are  bypassed 

No  Thrust 

Thrust  is  assumed  zero  and  mass  pro- 
perties are  constant 

No  CG 

Translation 

No  moment  translation  dxie  to  CG  off- 
set is  considered 

No  a , 6 
required 

The  aerodynamic  moments  are  not  a 
function  of  a or  6 

on  the  maximum  allowable  problem  dimensions  given  by  Table 
9.  The  two  test  cases  given  in  the  table,  discussed  in  de- 
tail in  Section  V,  provide  examples  of  the  algorithm  execu- 
tion times.  Execution  of  test  case  1 requires  approximately 
7 minutes  while  test  case  2 required  17.5  minutes. 


TABLE  9.  PROGRAM  SPECIFICATIONS 


PBOnRAM  spec:  FU' AT  I ON' 

MAXIHIIM 

ALbOHAnLF. 

TKST 

CASH 

3 

TEST 

CASE 

n 

Nijrrtji'r  of  dyiirimir  stntP.'; 

IB 

IB 

18 

Nunilmr  nf  et>tin»nlt*d  paTann'^tfrs 

.SO 

1 1 

:io 

Number  of  meaKiircmotitK  a1  e;iv.h  sample  time 

10 

R 

B 

Number  of  measuremnnt  sani.-'es 

6.B0 

100 

100 

Intecration  stepr  per  mr;i.,.irempn  t 

X 

1 

1 

'Maximum  allowable  dimensions  are  easily  modified  by  recom- 
piling the  algorithm  source  code. 
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Previous  experience  with  conversion  of  computer  al- 
gorithms similar  to  the  one  proposed  here  to  a CDC  6600  com- 
puter indicate  an  improvement  in  execution  time  of  at  least 
a factor  of  five.  Core  requirements  should  be  within  120K 
CDC  6600  words . 


SECTION  V 

ALGORITHM  VERIFICATION  WITH  SYNTHETIC  DATA 


INTRODU'CnOK 

The  alfioi'ithin  ver il  icai ion  results  presented  in  this 
section  are  based  on  the  processinp  of  compu i er-generated 
synthetic  rneasurerrient  s from  an  airframe  model  of  an  existing 
short-range  air-  t ..'-air  intercept  missile.  The  primary  pur- 
pose oi  fhos'-  scudies  i.s  to  demonstrate  (.he  validity  of  the 
extendv'd  ILilr'-oi  filter  as  an  identification  algorithni,  with 
a controlled  tost  case.  The  actual  state  and  parameter  esti~ 
marion  errors  are  compared  to  the  statistics  of  these  erioi s 
as  predicted  by  the  algorithti.  Consistency  between  the  actual 
error.s  and  the  predicted  statistics  is  an  indication  of  proper 
filler  perf ormanco , 

Divergence  of  the  estimation  errors  from  the  filter 
predicted  error  standard  derivati''ns  can  result  from  several 
sources,  stsch  as; 

• Inconsistency  between  the  measurement 

generator  model  and  the  filter  design 
model 

« Analytical  or  programming  errors 

• Nuiit  . leal  precision  (word  length) 

• Inaccurate  linearization. 

The  first  error  source  is  removed  from  consideration  in  the 
controlled  experiment  discussed  here  — certain  subprograms 
used  within  the  filter  algorithm  are  shared  by  the  measure- 
ment gener.ation  routine  insuring  model  consistency. 


48 


Because  of  the  complexity  and  generality  of  the  dy- 
namic state  equations,  the  major  source  of  error  has  been  with 
the  analytical  derivation  and  programming  of  the  analytical 
partial  derivative  calculations.  The  numerical  partial  de- 
rivative calculation  provides  an  iridependent  verification  of 
the  analytical  computations. 

Filter  inaccuracies  caused  by  numerical  precision  are 
frequently  encountered  in  large  scale  Kalman  filtering  prob- 
lems. Special  purpose  filLering  equations,  such  as  square 
root  filter  algorithms  (Reference  21,  are  often  implemented 
solely  for  the  purpose  of  avoiding  word  length  di f i icul ties . 
All  results  presented  in  this  report  were  performed  using  IBM 
360  single  precision  arithmetic.  Several  identification  fail- 
ures have  occurred  in  instances  with  no  process  noise  and  low 
measurement  uncertainty  — resulting  in  highly  accurate  param- 
eter and  state  estimation  accuracy.  These  failures  were  at- 
tributed to  precision  difficulties.  With  normal  measurement 
noise  and  process  noise  inputs  no  precision  difficulties  have 
arisen . 

Filter  divergence  resulting  from  violation  of  the 
linearity  requirement  is  a common  problem  in  the  application 
of  extended  Kalman  filters  (Reference  12).  This  problem  can 
sometimes  be  solved  by  using  accurate  initialization  methods, 
thf'  injection  of  process  noise,  or  resorting  to  more  complex 
filters  based  on  second  order  Taylor  series  expansions,  sucl: 
as  the  second-order  Gaussian  filter  (Reference  12).  F.  ctu- 
nately,  the  cases  investigated  here  have  d<..monsi  i aei’  .-eiii.irk- 
able  linearity  for  large  parameter  and  initi;.]  uncer- 

tainty. This  result  has  been  verified  through  i'-'  numerical 
partial  dei-ivative  computations  mentioned  previously. 


i 
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A distinction  between  two  types  of  models  must  be  made 
in  the  discussion  of  filter  evaluation  with  synthetic  data. 

One  model  is  used  for  generating  synthetic  measurements  while 
a second  model  is  used  within  the  filter.  For  the  majority 
of  the  discussion  in  this  section,  these  models  have  identical 
structure  — yet  differ  by  parameter  values,  state  initial  con- 
ditioii  values,  and  process  noise  inputs. 

The  short-range  air-to-air  interceptor  model  used  for 
generating  measurements  was  designed  to  exercise  the  various 
modeling  options  discussed  in  Sections  II  and  III.  Its  mass 
charact'  ri  cl  ics , thrust  model  ^.nd  aerodynamic  coefficients  are 
typical  of  existing  missile  configurations.  The  coefficient 
values  chosen  for  measurement  generation  were  selected  by  a 
Gaussian  random  number  generator  with  mean  values  equal  to 
representative  wind  tunnel  values  and  standard  deviations  of 
about  25  percent  of  the  mean.  The  resulting  aerodynamic  model 
corresponds  to  no  actual  missile  but  provides  a suitable  demon- 
stration case.  The  airframe  has  no  pitch  rate  or  acceleration 
feedback  for  stability  augmentation  and  thus  stable  missile 
flight  can  be  achieved  by  commanding  open  loop  aerodynamic 
control  surface  (canards  in  this  case)  deflections.  The  roll 
aerodynamics  were  not  modeled  — aerodynamic  roll  coefficients 
were  fixed  at  zero  — and  yet  considerable  roll  motion  was  in- 
duced from  the  roll  perturbation  moments.  The  linear  pitch, 
yav;  and  drag  aerodynamic  coefficients  were  assigned  values 
neai’  the  low  a ngi '’-of -attack  wind  tunnel  values  for  the  short- 
range  interceptor. 

The  data  shown  in  this  section  were  taken  during  the 


first  two  seconds  of  missile  flight.  The  missile  velocity  j 

variation  is  fron.  Mach  1.1  to  Mach  1.6  during  this  time. 

Time-varying  thrust  :-nd  mass  characteristics  are  modeled 

during  this  period.  ■ 
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A NOMINAL  TEST  CASE 


A test  case  is  discussed  in  detail  in  this  section 
in  order  to  demonstrate  the  algorithm  output.  Subsequent  sub- 
sections present  variations  to  this  case. 

The  aerodynamic  model  used  consists  of  11  basic  pitch, 
yaw,  and  drag  linear  aerodynamic  terms.  If  the  missile  dynamic 
response  can  be  restricted  to  within  a suitably  small  region 
about  some  known  reference  condition  during  the  period  of 
identification,  this  parameter  set  will  accurately  represent 
the  airframe  response  without  additional  nonlinear  terms. 

This  set  of  coefficients  is  widely  used  for  missile  perform- 
ance studies  and  autopilot  design  studies,  and  values  for 
these  terms  are  relatively  well  known  prior  to  fli,_,ht  test. 

We  assume  during  this  test  that  these  basic  aerodynamic  terms 
are  known  to  within  about  25  percent  (one  sigma).  These  11 
parameters,  the  values  used  for  filter  initialization,  and  the 
initial  parameter  uncertainty  assumed  within  the  filter  model 
are  given  in  Table  10.  Ihe  "true"  parameter  values  used  for 
measurement  generation  (also  shown  in  Table  10)  were  computed 
from  a Gaussian  random  number  generator  using  statistics  iden- 
tical to  the  filter  model.  The  resulting  aerodynamic  force 
and  moment  model  is  given  by  Table  11. 

The  test  flight  simulated  here  is  representative  of 
an  air  launched  flight  from  about  22,000  feet  altitude  at  Mach 
1.1.  The  mean  initial  dynamic  state  values  and  their  standard 
deviations  are  given  in  Table  12.  The  initial  state  uncertain- 
ties are  realistic  for  flight  tests  of  the  airframe  studied  here. 

The  aerodynamic  control  surface  deflections  for  the 
nominal  case  were  not  designed  specifically  for  the  purpose 
of  aei-ouynamic  parameT.er  estimation.  The  pitch  and  yaw  control 


TABLE  10.  FILTER  INITIALIZATION  AND  PERFORMANCE 


TAT-LE  11 


NOMINAL  TEwST  CASE  AERODYNAMIC  MODEL 


AEKODYNAMIC 

oirrpirr 

MATHEMATICAL  MODEL 

, 

X Body 

Force 

AOC 

y Brnly 

Force 

AOtCypB  + Cy^^Ar) 

•,  Body 

Force 

Monvent 

About 

X 

0 

Mompnt 

About 

* 

AI>OfC%i  a + Cii,  + Cu  0 

-a  'CM  • “q 

-'M 

Moment 

Abt.''  u 1 

7 

MF.>(Cf^pP.  + CNjqAr  + CN^r 

5^^ 

TABLE  12.  NOMINAL  INITIAL  DYNAMIC  STATE  STATISTICS 


STATE 

imiTS 

MEAN 

STANDARD 

DEVIATION 

Body  Velocity  (x) 

f t / see 

1160.0 

100  0 

Body  Velocity  (yl 

f t/see 

0.0 

20. 0 

Body  Velocity  (z) 

ft /sec 

0.0 

20.  0 

Inertial  Position  (x) 

ft 

0.0 

200.0 

Inertial  Position  ( y 1 

ft 

0.0 

200.0 

Inertial  Position  iz) 

ft 

-22100 

200.0 

Df'dy  Rotation  Rate  About  x 

rad/sec 

0.0 

0.20 

Body  Bn“;atlon  Rate  About  y 

rad /sec 

o.o 

0.20 

Body  Rotation  Rate  About  z 

rad /sec 

0.0 

0.20 

Roll  Euler  Angle  e 

rad 

0.0 

O-O-J 

Pitch  Euler  Angle  6 

rad 

0.0 

0. 04 

Heading  Euler  Angle  u 

rad 

0.0 

0.04 

1 


deflections  were  extracted  from  simulation  of  an  air-to-air 
target  engagement.  The  pitch  and  yaw  control  histories  are 
shown  in  Figures  5 and  6.  An  operationally  required  0.4- 
second  control  input  delay  is  indicated  on  these  figures. 

Six  process  noise  inputs  were  included  as  perturba- 
tion inputs  to  the  force  and  moment  equations  in  generating 
the  .synthetic  measurement  data.  Noise  was  processed  through 
low  pass  filters  to  generate  time-functions  representative  of 
the  appropriate  correlated  process  noise  inputs.  The  pertur- 
bation inputs  are  specified  by  their  noise  bandwidth  and  rms 
magnitude.  A bandwidth  of  10  rad/sec  (0.1-seccnd  correlation 
time)  was  selected  for  all  perturbation  inputs.  The  noise 
entering  the  acceleration  equations  (r^)  is  given  the  same 

rms  acceleration  level  as  the  acceleration  measurement  un- 

2 

certainty  — 2 ft/sec  (0,06  g).  The  moment  perturbations 
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Figure  5. 


Pitch  Control  for  Nominal  Test  Case 


TIME  (»«c) 


Figure  6.  Yaw  Control  for  Nominal  Test  Case 

are  given  rms  magnitudes  of  50  ft-lb  in  pitch  and  yaw,  and  10 
ft-lb  in  roll  — roughly  equivalent  to  a 0.3-degree  rms  noise 
level  on  control  deflection  measurement  or  a 5-ft/sec  wind 
gust  rms  level.  A typical  realization  of  the  perturbation 
input  is  given  by  Figure  7.  This  particular  sample  was  input 
as  the  body  x-axis  force  perturbation  for  the  nominal  case. 

Measurements  considered  for  the  nominal  test  case 
were  three  accelerometers  located  near  the  missile  nose,  three 
body  rate  gyros  and  the  a-f?  sensors.  White  noise  sequences 
were  used  for  the  measurement  error  models.  The  noise  levels 
(one  sigma)  are  given  in  Table  13.  The  measurement  sample 
rate  was  50  hertz  for  this  case. 

TABLE  13.  NOMINAL  MEASUREMENT  NOISE  LEVELS 


MEASUREMENTS 

UNITS 

S i tv'il) 

DEVI.-.?;  tn: 

All  Accelerometers 

ft/ sec^ 

T.O 

All  Rate  Gyros 

rad/ sec 

0.02 

a-y  Sensors 

rad 

0.  02 

Figure  7.  Typical  Perturbation  Force  Input 
For  the  Nominal  Test  Case 

The  estimates  of  the  aerodynamic  coefficients  used  to 
initialize  the  identification  algorithm  must  be  approximately 
correct  in  order  for  the  EKF  to  provide  useful  results.  The 
required  accuracy  of  the  a priori  model  cannot  be  generalized 
and  will  vary  from  case  to  case.  Nevertheless  any  obvious 
errors,  .such  as  control  input  sign  erroi  s , etc,  should  be  re- 
moved prior  to  attempting  parameter  identification.  Assist- 
ance in  initial  model  checkout  is  provided  by  simply  integra- 
ting the  a p’-jori  model  and  comparing  the  result  with  the 
flight  rif'-c-n;]  emo’it  s . The  comparsion  for  the  nominal  test 
case  discussed  abo- o is  presented  in  Figure  8.  The  trajec- 
tory mismatch  indicated  by  this  figure  results  from  initial 
condition  errors,  parameter  errors,  and  the  process  noise 
inputs . 


The  roll  rare  response  resulting  from  the  process 
noise  input  utilized  for  measurement  generation  can  be  seen 


from  Figure  8(d).  The  roll  angle  resulting  from  this  case 
was  approximately  560  degrees  after  2 seconds  --  realistic  for 
an  airframe  without  roll  stabilization.  The  predicted  roll 
rate  is  zero  --  indicating  no  a priori  information  about  the  roll 
response  is  assumed.  The  predicted  pitch  and  yaw  rotation  rate 
responses  exhibit  large  deviations  from  the  actual  as  a result 
of  process  noise  and  parameter  uncertainty.  This  rotation  rate 
sensitivity  to  parameter  errors  indicates  that  measurements  of 
rotation  rate  should  be  useful  for  parameter  identification. 

The  angl e-of-attack  and  sideslip  responses  resulting 
from  the  control  input  are  larger  than  those  from  a more  real- 
istic model  of  this  airframe  because  of  the  omission  of  the 
nonlinear  aerodynamic  effects.  Nevertheless  the  extreme  varia- 
tion in  angular  motion  offers  an  excellent  test  of  the  extended 
Kalman  filter  ability  to  identify  parameters.  Identification  in 
situations  having  different  angle-of-attack  and  sideslip  time- 
histories  is  discussed  in  the  subsections  to  follow. 

The  interpretation  of  parameter  identification  results 
and  the  verification  of  accurate  filter  performance  is  aided  by 
the  visualization  of  a variety  of  algorithm  outputs.  The 
following  variables  are  available  from  the  EKF  as  functions  of 
time . 


• State  estimate  and  state  estimation 
accuracy 

• Parameter  estimates  and  parameter 
estimation  accuracy 

• EKF  residual  process  and  its  one 
sigma  bounds. 


In  addition,  parameter-to-parameter  correlation  coeffi- 
cients are  useful  since  individual  parameters  may  be  highly 
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Figure  8. 
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unobservable  and  yet  highly  correlated  --  indicating  some 
linear  combination  of  parameters  nuiv  be  estimated,  Parameter- 
to-dynamic  state  correlations  are  also  available  and  nay  allow 
insight  into  additional  measurements  which  would  produce  better 
parameter  observability. 

When  synthetic  measurements  are  processed,  the  actual 
parameter  and  state  values  are  available  as  well  as  the  filter 
estimates.  Thus,  actual  state  and  parameter  estimation  errors 
may  be  computed  and  compared  to  the  predicted  one  sigma  error 
bounds  as  functions  of  time.  The  parameter  errors,  state  errors, 
residuals  and  all  respective  one  sigma  bounds  for  the  test  case 
discussed  above  are  contained  in  Appendi.x  A. 

The  results  in  Appendix  A must  be  interpreted  In  a 
statistical  sense  — recall  that  a normally  distributed  ran- 
dom quantity  is  within  one  sigma  G8.3  percent  of  the  time, 
within  two  sigma  95.5  percent  of  the  time,  and  within  three 
sigma  99.7  percent  of  the  time.  The  data  in  Appendix  A rep- 
resents only  one  sample  generated  by  an  initial  starting  seed 
to  a ratidom  number  generator.  A more  thorough  statistical 
analy.sis  requires  a l.arge  number  of  random  trials  and  the 
collective  interpretation  of  the  ensemble  results.  This 
monte  carlo  type  of  1 liter  evaluation  has  not  been  performed 
because  of  the  expense  of  computing  several  monte  carlo  sam- 
ples. Nevertheless  the  single  .sample  results  for  state  and 
parameter  estiniation  contained  In  Appendix  A are  a useful 
indication  of  acceptable  filter  psrfoimance . 

Ai  additional  performance  indicator  lies  in  the  white- 
ness of  the  filter  residuals.  Visual  verification  of  residual 
whiteness  and  correspondence  of  the  residual  magnitudes  with 
tneir  one  signia  bounds  shows  excc'llent  filter  performance  as 
seen  in  Figure  A-3  of  Appendix  A. 
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A useful  summary  of  the  parameter  estimation  perform- 
ance for  the  nominal  case  is  given  in  Table  5.  This  table 
includes  xrue  parameter  values  , final  filter  estimates,  final 
filter  .‘^xandard  deviations  and  a figure  of  merit  referred  to  as 
n^.  A value  of  greater  than  about  three  indicates  a highly  im- 
probable parameter  error  and  thus,  possibly  poor  filter  per- 
formance. The  value  of  is  computed  according  to: 


where 

a^  = estimate  of  i^  ' parameter 

t h 

a^  * actual  value  of  i parameter 

= filter  predicted  standard  deviation 

In  practical  applications  where  the  actual  model 
structure  and  parameter  values  are  unknov/n,  good  correlation 
between  the  predicted  measurements  from  the  filter  and  the 
actual  measurements  is  an  indication  of  acceptable  filter 
operation.  Yet  if  the  final  parameter  estimates  from  filter- 
ing are  inserted  into  an  otherwise  a priori  model,  a poor 
correlation  between  computed  and  actual  measurements  may 
result.  Small  levels  of  emitted  disturbance  inputs  (process 
noise)  or  erroneous  state  initial  conditions  may  have  a sig- 
nificant effect  when  integrated  uver  many  samples.  This 
phenomenon  is  demonstrated  using  the  nominal  test  case.  The 
computed  pitch  rate  agreement  with  the  measured  pitch  rate 
shown  by  Figure  9 is  not  significantly  improved  when  com- 
pared v/ith  Figure  8(e). 


Figure  G.  Computed  PitcJi  Rates  Using  Final  Parameter 
Estimates,  A Priori  Initial  Conditions  and 
no  Process  KoAse 

An  alto'cna  t )-ve  procedure  can  be  used  to  indicate  im- 
proved parameter  accuracy  when  synthetic  data  is  used.  In 
this  case,  state  initial  conditions  and  process  noise  inputs 
are  known  perfectly.  Thus , a priori  parameter  values  can  be 
used  with  perfect  state  initial  conditions  and  known  disturb- 
ance inputs  to  demonstrate  the  model  uncertainty  resulting 
from  initial  parameter  errors  alone.  Figure  10  shows  the 
pitch  rate  measurement  agreement  in  this  case.  A similar 
tra.iectory  produced  using  the  final  parameter  estimates  (with 
perfect  state  initial  conditions  and  process  noise)  resulting 
from  filtering  should  show  improved  agreement — indicat ing  model 
improvement.  This  is  demonstrated  in  Figure  11  where  excellent 
agreement  is  obtained  between  actual  and  predicted  pitch  rate 
measurement.  Although  not  included  here,  similar  excellent 
agreement  occurs  for  the  remaining  measurements. 
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Figure  10.  Computed  Pitch  Rate  Using  A Priori  Parameter 

Values,  Perfect  Initial  Conditions  and  Correct 
Process  Noise 
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Figure  11.  Computed  Pitch  Rato  Using  Final  Parameter 

Values,  Perfect  Initial  Condition.s.  and 
Correct  Process  Noise 


An  alternate  model  will  be  used  as  a reference  in 
certain  sensitivity  studies  to  be  shown  in  the  sections  to 
follow.  This  alternate  nominal  case  includes  an  erroneous  CG 
transl;tion  term  and  yet  the  resulting:  sensitivity  studies 
provide  useful  performance  data.  An  erroneous  sign  of  the  CG 
translation  term  resulted  in  a more  stable  airframe  and 
smclirr  angle-of-attack  and  sideslip  responses  to  the  control 
inputs  given  in  Figures  5 and  G.  The  anglc-oJ -attack  and 
sideslip  responses  for  this  case  are  shown  in  Figures  12 
and  13.  The  initialization  and  measurement  uncei'taint  ies 
are  iden''jcal  to  those  used  previouslt  . A summary  of  the 
pararccier  identification  performance  is  given  in  Table  14. 
Note  that  the  juirametor  estimation  accuracy  prc'dic.ted  by  the 
filter  is  worse  in  all  cases  than  for  the  previously  defined 
nominal  case  (Table  10).  This  is  a direct  result  of  the 
reduced  a and  B responses  for  this  case  --  also  note  that  the 
values  of  n in  most  cases  are  below  those  in  Table  10. 
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FLIGHT  SENSOR  STUDIES 

Adequate  flight  instrumentation  is  necessary  for  the 
extraction  of  useful  aerodynamic  models  from  missile  flight 
tests.  Three  studies  are  discussed  in  this  subsection  relating 
to  flight  sensors.  Studies  of  this  nature  £.re  useful  for  the 
design  of  flight  test  maneuvers  and  instrumentation  specifications. 

SENSOR  NOISE 

Sensor  errors  can  be  classified  into  two  categories 
for  the  purpose  of  error  modeling  — random  effects  with  no 
time- correl ation  called  "white"  noise,  and  errors  which  are 
correlated  in  time,  White  noise  is  easily  modeled  for  stat- 
istical filtering,  while  the  correlated  errors  require  a 
thorough  understanding  of  the  sensor  dynamics. 

Test  cases  for  four  noise  model  variations  from 
the  nominal  case  are  described  in  Table  15.  A summary 

TABLE  15.  NOISE  SENSITIVITY  STUDIES 


CASE 


DESCRIPTION 


1 


Two  times  nominal  noise  modeled  within  the  filter 
— two  times  nominal  noise  used  .for  measurement 
generat ion . 


2 


Two  times  nominal  noise  modeled  by  the  filter 
— nominal  noise  levels  used  for  measurement 
generation . 


3 Half  nominal  noise  modeled  in  the  filter — half 
nominal  noise  used  for  measurement  generation. 

4 Half  nominal  noise  modeled  in  the  filter — 
nominal  noise  used  for  measurement  generation. 


t 


? 
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of  the  filter  performance  for  the  first  two  cases  after 

2 seconds  of  filtering  is  given  in  Table  16  along  with  the 

nominal  filter  performance.  The  parameter  error  is  increased 

in  both  cases.  In  addition,  values  of  n are  reduced  from  \ 

^ i 

the  nominal  case,  indicating  better  predicted/actual  filter  1 

correlation.  Typical  residual  processes  from  these  two  cases  ' 

are  shown  in  Figures  14  and  15.  In  Figure  14,  where  the  filter 
model  is  correct,  64  percent  of  the  residual  points  are  included 
within  the  one  sigma  bound.  In  Figure  15,  96  percent  of  the 
residual  points  are  within  one  sigma.  This  result  demonstrates 
how  measurement  noise  statistics  can  be  identified  by  inspec- 
tion of  the  residual  process  and  the  predicted  residual  bound. 


I 
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TABLE  16.  PARAMETER  IDENTIFICATION  PERFORMANCE  WITH 

T Mr'nTr  A C’T'rx  a C‘TTr>T'»ir:»XTrn  TTrsTOT? 


PARAMETER 

NOMINAI, 

CASE  1 

INCREASED 
MEASUREMENT  NOISE. 
CORRECT  FILTER  MODEL 
o(n^) 

CASK  2 

NOMINAL 

MEASUREMENT  NOISE, 
INCORRECT  FILTER  MODE). 
o(n^) 

'=■'0 

0.0381 

(1 .6) 

0 . 0523 

(1.05) 

0.0541 

(0.99) 

0.574 

(1.6) 

0 899 

(0.81) 

0.929 

(0.67) 

0.00448 

(0.2) 

0 00586 

(0.05) 

0.00583 

(0.10) 

C 

0.692 

(1.5) 

1.07 

(0.78) 

1.11 

(0,67) 

q 

0.00601 

(1.1) 

0.00767 

(0.81) 

0,00766 

(0.56) 

Si 

a 

3.44 

(1.3) 

4.73 

(0.81) 

4.91 

( 0 . 66 ) 

Si  6 

q 

0 . 0843 

(1.07) 

0.1061 

(0.78) 

0. 1087 

( 0 . 63  ) 

q 

1032 

(0,010) 

1061. 

(0.04) 

1057 

(0.15) 

‘S: 

3.30 

(1.9) 

4.64 

(1.08) 

4.78 

(0.90) 

Si 

0 . 0608 

(1.3) 

0.0804 

(0.58) 

0.0819 

( 0 . .57  ) 

850 

(2.2) 

878.0 

(2.2) 

872.3 

(2.2) 

ASUREMENT  ERROR  (radi  YAW  RATE  MEASUREMENT  ERROR  (radi 


The  two  sen:.or  noise  cases  with  reduced  noise  levels 
are  summarized  in  Table  17,  and  compared  with  the  nominal 
estimation  performance.  For  these  cases,  0.8  second  of  data 
were  processed.  Parameter  accuracies  are  improved  in  both 
cases.  Note  that  the  values  of  are  larger  than  for  the 
nominal  ca.se. 


Numerical  difficulties  were  encountered  with  Case  4 
after  1 second  of  integration,  because  the  predicted  residual 
covariance  matrix  became  ill-conditioned.  Both  state  and 
parameter  estimates  show  signs  of  filter  divergence  prior  to 
this  point.  (Note  that  relatively  large  values  of  for 
this  case  at  0.8  second).  The  residual  process  for  the  angle- 
of-attack  measurement  from  Case  4 is  shown  in  Figure  16. 

Only  33  percent  of  the  residual  points  are  included  within 
the  one  sigma  bound. 


Several  characteristics  of  the  extended  Kalman  Filter 
are  demonstrated  by  the  noise  sensitivity  examples  discussed 
below : 


• Predicted  and  actual  estimation  error 
agreement  is  improved  with  increasing 
measurement  noise  levels. 

• Conservative  measurement  noise  levels 
have  little  effect  on  the  resulting 
filter  estimates  and  result  in  stable 
filter  performance  (Reference  12,  p.  244- 
251).  Correct  noise  levels  can  be 
inferred  from  the  resulting  residual 
process . 

• Low  measurement  noise  levels  can 

result  in  filter  divergence  and  numerical 
difficulties.  (This  could  be  corrected  by 
use  of  higher  precision  arithmetic). 


TABLE  17.  PARAMETER  IDENTIFICATION  PERFORMANCE  WITH 
REDUCED  MEASUREMENT  NOISE 


I'AHAMITKR 

9 

CASE  3 

REDCCED 

MEASUREMENT  NOISE. 
CORRECT  FILTER  MODEL 
o(n_^) 

CASE  4 

NOMINAL 

measurement  noise, 
INCORRECT  FILTER  MODEL 
o(n^) 

C 

■'"l. 

0. 110 

(0.27) 

0.0727 

’ .2) 

0.07.12 

(2.4? 

1.88 

( 0.1) 

1 . 32 

(1.6) 

1.310 

(2.8) 

0.00791 

(0.53) 

0.00780 

(0.7) 

0.00770 

(0.08) 

2.227 

(0.48) 

1.53 

(1.1) 

1.525 

(2.3) 

0.00872 

(0.54) 

0.00809 

(1.4) 

0.00823 

(0.55) 

s 

10.71 

(0.11) 

8.57 

(1.4) 

8.76 

(2.4) 

0.207 

( 0.29) 

0.1G9 

(1.6) 

0.176 

(2.4) 

2903 

( 0.38) 

2860 

(0.5) 

2867 

(0.35) 

7.82 

( 0.05) 

7.32 

(1.6) 

7.50 

(2.5) 

c 

1.83 

(0.80) 

0.155 

(1.1) 

0.160 

(1.63) 

r 

2598 

(0.33) 

2527 

(0.52) 

2542 

(0.50) 
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Figure  16.  Residuals  for  Case  4 


SENSOR  SELECTION 


Sensors  frequently  used  for  gathering  flight 
tesi  data  are  d.Lsctssed  in  Section  III.  The  generality  of 
the  extended  Kalman  Filter  facilitates  transition  betwe^en 
sensor  sets,  providing  the  appropriate  measurement  models  and 
partial  derivatives  are  included.  As  a demonstration  of  this 
sensor  selection  capability,  identification  results  from 
different  sensor  sets  are  presented  in  this  subsectioii. 

These  sensor  sets  are  given  in  Table  18. 

Case  5 is  the  same  as  the  nominal,  except  for  the 
longitudinal  accelerometer.  In  many  cases  missiles  have 
operational  requirements  for  only  two  accelerometers  mounted 
orthogonal  to  the  missile  centerline.  Thus,  flight  tests 


TABLE  IS.  SENSOR  SELECTION  STUDIES 


CASE 

DESCRIP'TION 

5 

Omit  longitudinal  accelerometer  from  nominal 

measurement  set 

6 

Accelerometer  measurements  only 

7 

Rate  gyro  measurements  only 

8 

Position  and  attitude  measurements  only 

using  only  operationally  required  sensors  would  not  include  a 
longitudinal  accelerometer.  The  primary  utility  of  this 
acceleromeTer  i.r  parameter  identification  is  in  the  deter- 
mination of  the  drag  related  aerodynamic  terms  (e.g.,  in  the 
linear  airframe  model). 

Estimation  performance  for  Case  5 is  shown  in  Table 
19.  The  predicted  drag  estimation  accuracy  is  degraded 
from  the  nominal  case;  however,  the  initial  drag  uncertainty 
of  0.4  was  improved  to  0.17  --  indicating  some  drag  informa- 
tion is  available  without  a longitudinal  accelerometer. 

Drag  information  is  apparently  inferred  from  know- 
ledge of  the  velocity  profile  — velocity  is  ''learned"  along 
with  the  dynan’ic  jiressure  through  the  known  relationship 
and  uncertaiiiLiGs  governing  the  lateral  accelerometer  mea- 
surements, dynamic  pressure  and  the  uncertain  parameters. 

The  initial  longitudinal  velocity  uncertainty  is  100  ft/sec; 
this  error  is  reduced  to  10,7  ft/sec  for  the  nominal  case, 
and  to  14.8  ft/sec  without  the  longitudinal  accelerometer. 

Accelerometers  and  rate  gyros  are  the  most  common 
sensors  found  on  operational  missiles.  In  Cases  6 and  7, 
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TABLE  19. 


PAPAMETER  ESTIMATION  ACCURACY  WITHOUT  A 
LONG I TUD INAL  ACCELEROMETER 


V'AKAw-.TilR 

NOMINAL 

0 { 11  „ ) 

CASE  5 

MO  L0N0!TUDINAL 
AOCELEROMETER 

0 ( .’> . ) 

0.''ia61 

a . P ) 

0.170  (0  IE' 

0.  sn 

(1.6) 

o.Gll  (3..1) 

S.r 

'0  CO-UK 

(0.20) 

0.00851  (O.l'i) 

C... 

0. 302 

(l.S) 

•3.782  (1.2.3) 

: _ 

0. OOCOl 

(1-1> 

C.OOeiU  (1.1) 

Si. 

3 . 4-i 

3.58  (1.1) 

0.08.13 

(1  OT) 

0.0862  (0.H3) 

3032.0 

(C.Ol) 

:o:'.5  (0.0.'!.' 

;) . 3u 

\0.01) 

3.86  (18) 

c . 

0,0608 

(1.3) 

0.0628  (1.1) 

'it 

S' 

800.0 

(2.2) 

b.S'J  (2  2) 

_»  

■ ... 

s"noor  yvti’  coi'.  i of  only  aocoi  oriiincifr.s  and  only  raif 
Kyro::i,  j'osi'Ocl  i vfly  . are  r.or..sidcrod . Tlu'  }>uramol.or 
mijon  ji'jri 'arni.M'.t’C  summury  Tor  those  two  cusos  is  shovf/i  in 
EG,  Lv.T  itiiatioi  aocu ru c i os;  and  values  oJ  for 
Case  C are  'shown  iil  0.8  tscvond.  ENtjc-me  fili.or  di verir.oi.ce 
Is  evjdoiC.  < as  ind.ic'itod  l>y  Itu  laijU  value.s  oi  ) \Cilch 
prevoiiled  yrocosslvi'j:  onta  beyL.-id  C 8 second.  No  lilier 
d i '.'ergon  tc  was  lound  wl.et-  only  rale  gyro.--.  :n  •.  us-'c  ami  the 
11  vnlueu  for  rlii.s  .'.KSe  indicate  fu-tei'  pei  1 'finance . 

Cu.‘.se  8 coinslslb  oJ  nieuou tnaicn t uJ  only  pJbilJ'.'ii  and 
u V.  l i ' ude  . T'.is  sol  iii'iy  result  IrNirii  jihoi  utlK.'odol  1 to  duLn,  ironi 
tt  iiiisblle  'Attliout  uaVio.ard  t elniii  I ry , or  j rorii  Iree  llly.lit  test.s 

i J 


TABLE  20.  PARAMETER  ESTIMATION  ACCURACY  WITH 
ACCELEROMETER  AND  RATE  GYROS  ONLY 


PARAMETER 

CASE  6 

ACCELEROMETERS 
ONLY;  T-0.8  sec 
o(n^) 

CASE  7 

RATE  GYROS 

ONLY;  T=2.0  sec 

0 (n  ) 

0.120 

(0.89) 

0.393 

(0.89) 

1.37 

(11.1) 

1.40 

(1  .R) 

*'  r 

0.00973 

(0.89) 

0.0098 

(0.55) 

C 

z 

a 

0.975 

(40.2) 

1.76 

(1.13) 

C 

0.00973 

(1.62) 

0 . 00996 

(0.26) 

Si 

a 

12.29 

(5.28) 

7.77 

(0.11) 

('<1 

0.228 

(^-8) 

0.192 

(0.49) 

S 

q 

2370 

(3.C) 

529 

(2.0) 

r 

R « O 

r 1 o n \ 

V ^ . » vy  / 

*7  rvi.' 

-Ne, 

« » 

V.  L'  . %JKJ  ’ 

Si 

r 

0.3  93 

(1.6) 

0.;i23 

(0.31 ) 

S: 

2780 

(0.73) 

5}  3 

tnv*  Mf^«w 

( 1 . 0.3  ) j 

oJ  srnu  1 1 -eculc  modolp.  flit-ti  j-i'O.ircl  1 Iv  test 

rji?;Les  (Udcrence  IS).  The  j>03’odyi;ari» J c;  lorcc*  and  niouK.Tl  co- 
ol ficicTits  (julvr  directly  Into  the  dynamic  oqui.  I i .lor 
l.iricar  uud  niiyului  Acci  Ic: 't;  loii---  tlius,  hccoI  era  l i on  Hiif.',  j;r)- 
Kulur  rute  provide  niori'  dii’ect  tnlorrnal.i.  OTi  'il, '.lU*  "‘.■I'l.t!  .. '.Mt'O I'^y '• 
nunii^’  1 ers  than  dt)  V)o;iii<,ion  iind  m.lili.iC(  n',ea/j’,;jv’W  vTi'i.j , 
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Nevertheless,  highly  accurate  position  and  attitude  measure- 
ments have  been  used  with  success  to  extract  aerodynamic  data 
in  ballistic  test  ranges. 

Three  missile  position  measurements  and  three  Euler 
angle  attitude  measurements  were  assumed  with  measuremen c 
noise  le^'cls  of  5 feet  and  0.0''  radian,  respectively,  at  & 
data  rate  of  50  hert".  Those  accuracies  are  representative 
of  phoiotheodol ito  measurements.  The  resulting  parameter  est 
mat  ion  accuracies  are  compared  v;ith  the  alternate  nominal  re- 
sults in  Table  21.  The  filter  performsuce  in  this  case  cor- 
related well  wiili  jm'edicted  peri ormarce , hovev'c*)’,  the  result- 
ing parameter  accuracies  were  degraded  from  th#?  noiyiiPRl  mea- 
surement set.  It  is  interesting  to  note,  in  coa.paring  Ccses 
7 and  y,  that  more  aerodynamic  information  is  available  from 
i-ato  gyros  alone  than  from  position  and  attitude  data,  for 
tpe  measuron-.'-nt  uocuruc  jes  lr,\  ost  igateU. 

A signl  t leant  isisue  witii  regarci  to  sensor  s<f  lection 
ct-ncerns  th»>  iiecs'sst  ty  of  the  a-n  seusi'l' . Tl.is  issue  will  he 
addvossed  in  ntvrv-  detail  later  i;  report.  The  alternate* 

ru.'imnn,;  »va;.;  repro'.:ca.i>: v.  using  only  acet  Itromeicr!-.  Riid 

rate  gyros.  No  difference  was  noticed  IP  tr»e  r<*f.ul i i yaram 
eier  i done  i f i a t ion  pr  r f ormunce  , indicating  little  additional 
i »i X orma I iosi  is  contributed  by  tii'-  o-h  sense.-;  . 

•''he’bO/!  «f0'0bLl Kd 

£t.  vernl  stuoics  were  perior.'i '''d  to  iin  cst  igi.t  v t h-. 
t-i  bybtemutJc  sens*,)/  errors  oc  v-y.r;v’i'.e; -.■r  iden t i f 1 cu- 
ti'.;-i  'J'l;- •■•■u;  oy;>i  wl'iUt  Jo  4 rj'on,  roelvu'b-  ),  (.  rulii<*  t c r 

and  bC.Fj'v;  Ijo  to)  crior,  and  v-.-.;.o  biu;,  ki  O oouie  l actor  e}i>w'r. 
V‘!'/i',h  :t.ie  inodc.lid  us  t':-, sea '.nod  J »i  j s><.  111.  The  Ciiues  pi  i- 
u here  -.lU.  'Mr*-  id.\:»it  , vai  meuauv cf^^eiC,  i;.  juh.  1 a und  IJltei 


TABLE  21. 


PARAMETER  IDENTIFICATIOK  PERFORMANCE  WITH 
POSITION  AND  AT'i'ITuDE  MEASURLMENTB 


r\vUcl?  so  ir.iit  wh<;n  l'li‘6  yrroTfi,  fur  ev.a.nj  le,  urc;  iiiclud<,*d 
In  i.bie  nioi-.-^urcrttnl ; Hioy  ar<'  ul.so  n'.odt-lG.I  and  oKt inia/’.c-cl 
\J:{.  3' j,  i \ . 

TL';  cl>nR»iic  kind  r-i^i' l 1 ''^v ici  1 tnocW.  i.  ui>f;d  in  y.onorai, in}/, 
il'.i'  k'.l -.03  nr-.tv  iio:iii)ih  J.  vii\jncv.ui-y  war  ui.ll1^t;’J  io  :An:  illLcr/> 
cliscg.wv-od  in  itilci  uubsvoi  i '--i! . Tho  bnnd  rn(,4Mt j.'oU  /or 

».;orior:i'.  IviS,  w)  '..  x luid'./in  (»t./;ya«'..’)v'xvv  \>rt<.‘u  ar  nioaw'jr/'tnoia  y.nntjx  v l ioji 
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vas  different  from  the  seed  used  to  produce  the  results 
foi'  the  alternate  nominal  case  shown  in  Table  14.  This 
results  in  different  perturoation  input  time  histories, 
parameter  initializations,  and  meacuvement  errors.  The 
alternate  nominal  parameter  estimation  results  wore  recom- 
puted (without  systematic  errors)  with  the  new  seed  ior 
comparison. 

The  bias  and  scale  factor  levels  studied  are  shown 
in  Table  22.  Each  type  of  error  was  studied  individually, 
lor  a total  of  four  expunples;  each  of  the  foui‘  cases  included 
the  three  additional  parameters  associated  with  the  corre- 
sponding .systematic  errors  — x,  y,  and  z axis  biases,  for 
■example.  The  predicted  p.-^rameter  identification  accuracy 
for  the  four  cases  is  shown  in  Ta)  .t.e  23,  together  with  the 
•alt eriiv.le  uoininal  case  having  no  tystematic  errors. 


SYSTEMATIC  ERRORS  SENSITIVITY  STUDIES 


•iiM'rHWWMc 

CASE 

DLPCRTPTIO’W  j 

NO. 

INihOR  TYPE 

^3TANDARD  DF.VJ  AT.TON  ] 

9 

Acc'-.leromet  or  bias 

2.0  I t / sec 

10 

Accelcroivetcr  .scale  i actor 
or  ro  r 

0 . 02 

I 1 

Kate  ijyio  bias 

0,02  !■  ad,  ' -c 

I 2 

31?.  t c (/  y r <j  c v.  1 o f a (.;  t o r 

error 

o.oi; 

TABLE  23.  THE  EFFECTS  OF  SYSTEMATIC  MEASUREMENT  ERROR 
ON  PARAMETER  IDENTIFICATION 


PRrnicTF.D  parameter  .standard  deviations* 


PARAMETER 

ALTERNATE .NOMINAL 
K.)  SYSTEMATIC 
ERP.On 

0A:;i:  9 

A'-’CELEROStirTEIi 
nlAS  EliRORS 

CASE  1 0 
ACCELT-ROMETE}! 
SCALE.  FACTOR 
..JiERUES,  ™ 

CASE  11 
. GYRO  BIAS 
ERRORS  1 

CASE  12 
GYRO  SCALE 

1 ACTOR  ERRORS 

i).0Z3l 

0.0387 

0.0370 

0.02-17 

0.0250 

0.  REP 

0.91V 

(1.87-1 

0,090 

1.021 

r 

.0.  OOHKG 

0.0067? 

0.00/45 

0.00572 

0.006P1 

f 

='u 

1.1  G9 

. • 4 

.1 . R09 

1.190 

1.0-10 

1 . 327 
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i “bj 

0.00028  j 
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O.OOC38 

0.00710 

0.00OR2 

o.o'invs 
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6 . 90 

0.  05)2  4 

0.0910 

[ 
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conti'ol  surface  orientation  with  respect  to  the  relative 
wind.  The  aerodynamic  coefficients  relate  the  reJative  wjud 
to  the  forces  and  moments.  To  estimate  the  aerodynamic  co- 
efficients, estimates  o."  force  and  moments  and  estimates  of 
the  relative  wind  must  be  inferred.  Given  the  aerodyncUtic 
coefficients,  the  relative  wind  information  (a  and  C)  can  be 
detei'mined  from  acceleration  and  rate  gyro  ineasureuient  and 
known  control  surface  del  lections , Also,  given  good  esti- 
mates of  initial  missile  inertial  velocity  and  a perfect  mea- 
sure of  atmospheric  v/ind  the  relative  wind  velocity  can  be 
estinsateci  via  integrs  ti  on  oX  the  acceleration  over  short  time 
intervals.  With  uncertain  wind  and  initial  velocity  errors 
as  well  as  parameter  unce.rtainty , it  i.s  not  clear  whether 
0 and  6 can  be  estimated  with  sufficient  accuracy  to  improve 
the  estimate.s  of  the  aerodynamic  parameter  . Thus,  it  might 
t-  argued  that  c/.  and  £ sensors  a.re  required  for  p.Ava.meter 
icient i .f icat ion  when  atmospheric  wind  uncoctaint .los  are  pr^,'sont. 
A brief  study  nddresoing  this  issue  is  discussed  here.  The 
test  case.s  are  listed  in  Table  24, 


I J.  A Oi  t:  1 uQ  J.  li  t:,  c?,  ^-n:;  ; m 1 j .;>  v 4.  v»,».a,v*  j. 

included  in  the  hKF  general.i.tod  airframe  model  . The  steady 
compoiient  oi  wind  is  moueled  by  a horizontal  vector  with 
magnitude  and  heading  that  arc  functions  of  altitude.  The 
inclusior-  of  a steady  wind  correctly  riiodelcd  within  the  filter 
dynatDiC.s  was  found  to  have  litt.le  c3  feet  on  the  j-'Cramoter 
ident  i f i ca\.  J Oh  perlormunoe  (a  slight  Identification  irr.prov*,'- 
nient  was  noticed,  tipparcnlly  fj’oni  the  additional  a and  v’^riu- 
tJo)i  result  jng  from  tht;  wind.  A nnrjzonta’.  v ; r.d  valv-  ol  about 
110  Xt/ucc*  oi  Icnted  45  dogicvfS  from  t he  in  i K' I mjv-slle  head- 
n.Mg  va,s  assumed. 


’♦hO  jJcrctiM.ile  Wind  ut.  20l<  leeL  alU’urde  (Kviereeve  14} 
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TABLE  24.  WIND  STUDIES 


CAii;- 

DESCRIPTION 

13 

Nominal  case  with  deterministic  wind  added 

Uuc-eri. aia  wind  added  to  Case  13 

15 

Ca.oe  14  with  poor  a-B  sensors 

1 

16 

Case  15  with  poor  initial  attitude  and  velocity  j 

va  iue.s 

17 

Ca.te  16  with  accurate  u-S  sensors 

S 


In  addition  to  tlie  cleterininxst ic  v/xnd,  a random  (but 
constant)  wind  biaa  v/as  niodeled  and  estimated  (Case  14), 

Random  v/ind  blase-;  o.i'  up  to  50  ft/sec  in  the  three  inertial 
ax'GS  were  modeled.  With  nominal  measurement.s  and  initial 
conditions,  little  change  v/as  observed  in  the  parameter  iden- 
tiiication  accuracy  and  the  wind  uncertainty  vras  improved  J.roro 
50  It /sec  to  about  SC  it/sec  in  each  oi  the  three  inertia] 
direction-s  - The  nominal  c and  g sensors  were  included  in 
these  cases. 

The  el  loots  ol  the  a aiid  0 sensors  on  the  idontili- 
cation  accuracy  with  larp.e  deterministic  wind  and  large  wind 
uncertainty  was  determined  by  varying  the  measuren cat  noise 
level  on  t/iO  u and  fc-  sen.sor:>  (Case  15).  The  nornin?.!  a and  1-. 
sensor  nois^e  Ic'e.l  ol  Q.02  rad.>.R.n  was  increased  to  0.2  radian, 
representing  es.seniially  no  itjlormatlon  provided  by  the'  a and 
t sensors.  Tor  this  cast,  there  was  no  noticubJ.e  Ut-'creuse  in 
the  paramoicr  Ident i 1 1 cation  accuracy.  boiiic  degradaLion  in 
the  ability  to  estimate  wind  v/ao  noticed,  with  wind  uncertain- 
ties ol  ab(.)uL  40  it. /sec. 


Some  additional  studies  were  conducted  to  isolate 


the  relationship  between  initial  state  uncertainties  and 
parameter  identification  in  the  presence  of  wind.  Note  that 
if  good  initial  estimate-s  of  the  lateral  components  of  mis- 
sile inertial  velocity  are  available  and  if  good  initial 
estimates  of  missile  attitudes  are  known,  good  initial  esti- 
mates of  a and  P can  be  computed.  By  propagating  these  ini- 
tial estimates  through  the  airframe  dynamics,  a and  6 esti- 
mates may  be  inferred.  For  Case  16,  the  lateral  velocity  com- 
ponent uncertainties  were  increased  from  the  nominal  value  of 
20  to  80  ft/sec  while  the  initial  missile  attitude  uncertain- 
ties were  increased  from  nominal  values  of  0.04  to  0.2  radian. 
With  the  large  a and  6 sensor  rms  noise  level  (0.2  radian), 
this  case  exhiblued  filter  divergence.  Figure  17  .shows  the 
resulting  z-axis  wind  estimation  errors  as  a function  of  time. 
Case  17  was  performed  with  verj^  low  u and  E sensor  noise  levels 
(0.002  radian);  the  filter  performance  was  near  nominal  with 
the  z-axis  wind  estimation  error  shown  by  Figure  18.  Thus,  to 
avoid  divergence  in  the  presence  oi  large  state  initial  con- 
dition errors  and  uncertain  wind,  accurate  a and  E measure- 
ments may  Vjo  required. 


ESTIMATION  OF  NONLINEAR  AERODYNAMIC  EFFECTS 

The  picvious  sections  considered  only  the  so  called 
"linear"  aerodynamic  effects,  however  the  resulting  dyuamit- 
model  i.s  highly  nonlinear  in  the  augmented  system  .state  vector. 
The  inclusion  ot  additional  "nonlinear  aerodynan.  lc  elJects" 
is  conceptually  no  more  difficult  to  handle  in  the  paiarneter 
identil  ication  aigoi’ithm.s  than  the  pievious  exumi  ies.  Never- 
theless, in  orde*r  to  demon.struie  the  caj  abJ  1 ities  of  identi  ly- 
ing a large  nurtiber  of  parameters,  a variety  t'f  nonlineur  and 
crouti  coupling  aerod/riwnic  parameters  ate  identified  in  this 
section . 
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Figure  17.  z-Axis  Wind  Estimation  Error  With 
Large  Initial  State  Uncertainty 
and  Large  u-B  Sensor  Errors 
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Figure  18.  z-Axib  Wind  Ebtimutl'm  Error  With 
Large  luiliuJ  Stale  Uncertainty 
and  Snail  a-f  Sensor  I'rrtirs 


One  strategy  for  applying  system  identification  to 

* 

an  airframe  might  include  the  following  steps.  First,  attempt 
to  fit  the  observed  data  using  only  the  linear  aerodynamic 
terms.  Observe  the  resulting  residual  process,  parameter 
estimates  and  the  respective  one  sigma  bounds  and  evaluate 
the  filter  performance  from  the  standpoint  of  consistency. 

In  many  cases  a linear  aerodynamic  model  should  provide  rea- 
sonable correlation.  Next  isolate  potential  nonlinear  aero- 
dynairdc  effects  likely  required  to  improve  the  model.  These 
additional  nonlinear  terms  may  be  provided  by  the  wind  tunnel 
data  or  through  insight  gained  by  inspection  of  the  residual 
process.  For  example,  if  a large  descrepancy  occurs  in  the 
lateral  acceleration  residual  at  a time  when  a large  roll 
rate  is  observed,  a roll  coupling  term  might  be  added.  This 
"intuitive  approach"  to  model  structure  selection  relies 
heavily  on  the  experience  of  the  user  and  requires  a highly 
flexible  identification  algorithm  so  that  aerodynamic  co- 
efficients can  be  removed  or  inserted  easily. 

Once  a list  of  potential  additional  parameters  has 
been  selected  to  improve  the  linear  model  results,  parameter 
uncertainty  standard  deviations  must  be  selected  to  initializc- 
the  additional  parameters.  In  some  cases  an  uncertainty 
level  might  be  ascertained  from  the  wind  tunnel  data.  Alter- 
natively, uncertainties  might  be  computed  by  considering  the 
contribution  of  the  additional  aerodynamic  terms  to  the  total 
aerodynaijiic  force  and  moment  coefficients.  For  example,  if 
the  average  angle--of-attack  and  sideslip  over  a segment 
selected  foi’  identification  is  observed  to  be  a and  B,  the 
average  contribution  to  from  a is  about  a.  The  uncer- 
tainly in  a term  such  as  might  be  selected  to  produce 

an  uncertainty  jn  C..  of  about  10  percent  of  the  known  i>riiiiary 
contributor,  . Thus,  the  standard  deviation,  , in 
IS  selected  as, 


I 


3 


0.10(Cj^  a)/(a  B) 


A St  was  conducted  where  25  additional  terms  were 
included  in  the  aerodynamic  model  used  in  the  nominal  cases, 
resulting  in  a total  of  36  aerodynamic  parameters.  These 
parameters,  their  mean  values,  standard  deviations  and  true 
values  are  given  in  Table  25.  The  standard  deviations  were 
selected  so  that  each  new  term  contributed  on  the  order  of 
5 percent  to  the  aerodynamic  coefficients. 


Different  pitch  and  yaw  control  inputs  were  used 
for  this  test  case.  These  inputs  together  with  the  resulting 
anglo-of-attack  and  sideslip  are  shown  in  Figure.^  19  and  20. 
Only  one  second  of  data  was  processed,  but  the  data  rate  was 
■>ncreasecl  to  100  hertz.  The  large  angle-of-ar tach  and  side- 
slip variations  in  this  trajectory  should  provide  good  ob- 
servability of  the  nonlinear  cerms. 


An  additional  test  case  with  36  parameters  was  per- 
formed where  the  control  inpuL.s  of  Figures  19  and  20  were 
reduced  by  half.  The  resulting  anglc-ot-attack  and^.si'd'e.slip 
are  also  indicated  in  Figures  19  and  20.  It  isi-iaterc  st  :i ng 
to  note  that  the  addition  of  the  noniaruM':!:'"  unci  ccMipljag  aero- 
dynamic effects  are  not  obviou^f^rOfi  casual  observation  of 
Figures  19  and  20;  the  lOitcb”  and  yaw  plane.s  do  not  appear 
coupled  and  the  a,  and  B responses  to  contio.'.  inputs  could 
result  from  a linear  aerodynamic  model.  The  idcntil ication 
o.f  subtle  aerodynamic  eflects  which  cont7lbut.e  little  to  the 
basic  linear  airfryine  .response  provides  an  excellent  demon- 
stration of  the  capabnities  cf  the  postllicht  data  pf^fos.sing 
algorithm  discussed  here. 
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TABLE  25.  PARAMETER  IDENTIFICATION  SUMMARY 
FOR  36  PARAMETER  CASE 
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Figure  19.  Angle-of-Atiack  Profile  fo: 

the  36  Parameter  Cases 
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The  predicted  parameter  accuracy  as  well  as  the  values 

of  n are  shown  in  Table  25  for  both  control  input  levels.  The 

o 

parameter  estimation  errors  and  the  predicted  one  sigma  bounds 
are  shown  in  Appendix  A,  along  with  the  residual  processes. 

All  consistency  tests  indicate  the  extended  Kalman  filter  to 
be  performing  correctly,  the  residuals  appear  white  and  re- 
sidual magnitudes  are  predicted  accurately  and  the  state  and 
parameter  errors  behave  according  to  the  predicted  bounds. 

The  values  of  n indicated  in  Table  25  indicate  excellent 
o 

filter  performance. 

With  the  larger  control  input  and  the  resulting  larger 
a and  B,  accurate  estimates  of  the  third-order  terms  are  ob- 
tained rapidly.  Note  from  Figure  A-4  that  the  primary  param- 
eter accuracy  improvement  comes  prior  to  about  0.4  second  -- 
after  about  0.5  second  little  additional  improvement  results. 
Yet  the  majority  of  the  large  a and  6 variations  occur  after 
0.5  second.  The  drag-related  terms  and  the  Mach  number  varia- 
tion effects  are  less  observable  than  the  pitch  and  yaw  dy- 
namic terms.  The  drag  terms  produce  relatively  slight  changes 
in  longitudinal  acceleration  so  that  difficulty  in  identifying 
these  effects  is  anticipated  over  short  data  segments. 

The  estimation  accuracy  achieved  with  reduced  control 
magnitude  was  significantly  reduced  for  ten  third-order  non- 
linear terms,  with  the  average  degradation  being  a factor  of 
about,  seven.  Estimation  error  degradation  for  the  remaining 
parameters  was  not  significant.  In  fact,  the  five  Mach  mu  .-er 
variation  terms  were  improved  with  less  a ant!  f-  var-.f.  ion. 


ADDITIONAL  PERFORMANCE  STUDIES 

A variety  of  additional  studies  have  been  conducted 
during  the  course  of  algorithm  development.  In  particular, 
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parameter  idont ii ica t ion  sensitivity  to  data  rale,  state 
initial  condition  magnitude,  and  process  noise  are  discussed 
in  this  section  (Table  26).  These  cases  are  variations  of 
the  alternate  nominal  trajectory.  The  predicted  identifica- 
tion accuracies  for  these  three  demonstration  cases  arc  given 
in  Table-  27. 


TABLb  26. 


ADDITIONAL  PERFORM.ANCE  STUDIES 


CASE 

NO. 

DESCRIPTION* 

1 

i ' 

Jv8 

Ilalt’e  the  nominal  data  rate  to  25  Hz 

1 ^ 

19 

Double  all  initial  dynamic  stave 

1 i 

uncer  k 'cl  j.  t)  i i 

n 

1 1 

20 

Set  the  process  noise  magnitude  equal 

i ! 

tC'  zero 

♦The  measurement  and  filter  model  structures  are 
identical  on  all  cases. 


Case  18  indicates  that  data  rate  reduction  increased 
the  paranieter  estimate  standard  deviations,  but  not  signifi- 
cantly. Doubling  the  sample  rate  (data  not  shown)  did  not 
change  the  estimation  accuracy  — indicating  50  hert:-:  to  be 
near  the  optimal  data  rate  for  this  nominal  case. 

Doubling  the  dynamic  state  initial  condition  uncer- 
tainty (Ca.se  19)  also  increases  the  parameter  estimation 
error,  but  not  s , r.cif  icantly . This  demonstrates  an  insen- 
sitivity of  the  estimates  to  state  initial  conditions,  as 
well  as  showing  that  large  initial  condition  uncertainty  can 
be  tolerated  without  indications  of  filter  divergence. 
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TABLE  27. 


ADDITIONAL  FILTER  PERFORMANCE  STUDIES 


iMi'i.M.  I’AHAMiTrn 
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24  7.0 

Ca.se  20  shows  i den  t.  i f .i  cat  ion  accuruc\-  with  no  pro- 
cess noise  input.  In  soiiie  cases  significant  identification 
improvement  results,  particularity  Cj^j^  aticl  . The  process 
noise  adds  directly  to  the  force  ard  moment  eciuations  which 
include  the  aerodynamic  parameters.  Ab  more  noise  is  added, 
the  deterministic  structure  of  the  monel  has  relatively  less 
effect  so  that  the  parameters  become  less  observable. 

An  additional  test  case  was  performed  in  order  to 
demonstrate  the  capability  to  fix  certain  well  known  param- 
eters while  estimating  other  terms.  The  less  visible  pa- 
rameters from  the  nominal  case  are  and  . 

If  the  remaining  7 coefficients  are  fixed  at  the  true  values 
those  less  observable  parameters  should  become  more  visible. 
The  improved  identification  accuracy  for  this  case  is  .-hown 
in  Table  28. 
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IprNTIFICATIGN  ACCLHArV  WITH 
y-'ARTIAL  T'ARAMETi-R  SET 


T-iSLE 


2S 


PAK/iVi  Ti 

IMTI/.I  I'NCMiTAIKTY 

4 I'ARAMETEPS 

^'v 

I ‘ ‘ 

O.OJ 

0. 00448 

0.00288 

■jij 

0.01 
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0.00432 

"V 

;U)0U  . (• 

lO.V  .0 

1020.0 

'■■'k 

aocio . 0 

Boo . 0 

79r> . 0 

1 

~ I.,,,  « ,.ll.  .. 

SYN'TdETIC  DAiA  PROCESSING  F'.'R  A CASE  EXHIBITING  LARGE  MACH 
NUMBER  VARIATIONS 

The  aitjoi-ithm  veriT icax jor?  exercise  presented  in  the 
previous  sections  uiilized  sisnulated  measurements  produced 
with  subroutines  sha  -ed  by  the  filter.  For  these  cases  the 
.structui'c-  of  the  dyu  imics  used  for  measurement  generation  was 
known  to  be  id.^nticai  to  the  strucv.ure  assumea  within  tne  fil- 
tering algorithm. 

An  addlitonaJ  study  was  to  be  performed  during  this 
GXfori:  in  order  to  e.xercise  the  genera]  purpose  airframe 
r.iociel  and  the  fi]tt;ring  algorithm  against  synthetic  measure- 
m'-'nt  data  j'rotjded  ny  the  Air  Force.  The  data  were  to  be 
geiioriitcd  a niissile  exhibiting  large  variations  in  Mach 
n’-imber . 


Synthetic  data  were  providea  by  AFATL  that  were 
generated  from  a six-degrees-of-freedom  simulation  developed 

♦Air  Force  Armainent  Laboratory,  Eglin  Air  Force  Bas^',  Florida. 
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indciiondEMit  ly  from  The  moilol  discussed  in  Section  H.  This 
data  rep.’.'esents  a G.G-aeoonds  trajectory  from  a high  lift, 
bank-to-turn  air  frame  for  a short-range  air-to-air  inter- 
cc'i’tor.  The  trajectory  contains  a 2.6-secoiids  boost  followed 
by  four  seconds  of  coast  flight  during  which  various  pitch 
and  yaw  maneuv<>rs  are  executed.  The  Mach  number,  angle-of- 
atiack,  and  sideslip  profiles  are  given  on  Figures  21,  22, 
and  23.  Other  tr.ajectory  variables  — accelerations,  body 
rates,  attitudes,  etc  --  were  provided  by  AFATL  to  represent 


noi se-f roe  measurements . 


The  equation.^  used  lor  modeling  the  aerodynamic  forces 
and  moments  were  provided  by  AFATL  and  are  shown  in  Table  29. 
Aerodynamic  parameter  values  corrupted  by  "approximately  10 
percent  error"  were  provided  by  AFATL  in  the  form  of  Table  30. 
The  objective  of  tins  exercise  was  to  identify  an  aerodynamic 
model  which  accur.ately  reproduced  the  synthetic  measurements 
during  several  segments  of  the  trajectory. 


i 


Time  (utri 

Figure  21.  Mach  Number  Profile  for  AFATL  Interceptor 
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TAULL  29.  AERODYNAMIC  MODEL  LOR  AFATL  INTERCEPTOR 


O.CVC'f'  I O.OOM 
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HOnVNAMir  PARAMETEK 
Tl,  INTERCEPTOR 
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The  significant  characteristic  of  this  data  from  the 
standpoint  of  parameter  identification  is  the  large  variation 
of  Mach  number  during  boost  and  coast  as  well  as  the  large 
variation  of  the  aerodynamic  coefficients  with  Mach  number, 
as  seen  in  Table  30.  The  percent  rate  of  change  in  param- 
eter magnitude  (PC)  at  about  Mach  2 is  also  indicated  on 
Table  30.  The  quantity  PC  is  computed  according  to, 


PC  = i — / C *100 

3M  at  / ^i 


where 

"t  h 

C^  = i ^ aerodynamic  coefficient 
M = Mach  number 
atl/dt  = 0.723  per  second 

Note  that  for  the  majority  of  parameters, PC  is  greater  than 
30  percent  per  second  for  this  trajectory. 

The  aerodynamic  parameterization  presented  in  Section  II 
assumes  only  slight  variation  of  the  aerodynamic  forces  and 
moments  with  Mach  number  during  any  flight  segment  selected  for 
identification.  The  allowable  Mach  variation  can  be  expressed 
more  quantatively  by  considering  the  time  required  for  accurate 
parameter  estimation  versus  the  time  variation  of  the  param- 
eters with  mach  number.  For  example  if  the  parameter  uncer- 
tainty at  Mach  2 were  10  percent  and  the  parameter  rate  of 
change  with  Mach  number  is  about  30  percent  per  second  — then 
the  allowable  time- interval  for  accurate  identification  should 
be  small  with  respect  to  1/3  second. 


It  becomes  apparent  that  the  aerodynamic  parameteriza- 
tion of  Section  II  is  inadequate  for  the  large  Mach  variation 
in  the  synthetic  data  provided  by  AFATL.  Several  alternative 
methods  are  proposed  to  allow  aerodynamic  modeling  from  flight 
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data  in  cases  where  all  parameters  may  have  large  but  line.r 
time  variations  over  the  selected  segment. 


• Estimate  the  parameters  as  though  they 
are  constant  — assume  the  resulting 
parameter  estimate  to  be  an  average 
value;  i.e.,  associate  the  final  esti- 
mate with  the  mach  number  at  the  center 
of  the  .segment. 

• Add  parameters  to  model  the  linear 
variation  of  each  coefficient  with  mach 
number;  this  will  double  the  number  of 
parameters  to  be  estimated. 

• Model  the  aerodynamic  coefficients  as  time 
varying  random  processes  --  i.c.  dynamic 
states  rather  than  parameters.  This  does 
not  increase  the  dimensions  of  the  required 
EKF.  Reference  15  provides  an  example  of 
identification  of  time-varying  parameters 
using  this  method. 

o Treat  the  aerodynamic  model  presented  in 
Section  1 1 as  a perturbation  model  which 
must  be  added  to  a more  thorough  a priori 
model  containing  the  primary  Mach  number 
variation  affects.  The  model  used  for 
identification  would  thus  include  the 
wind  tunnel  model  used  for  prefiight 
studies  and  would  generate  identical  sim- 
ulation results  if  all  perturbation  aero- 
dynamic coefficients  were  set  to  zero. 

The  perturbation  model  coefficients  remain 
time  invariant  and  yet  should  remain  valid 
over  larger  flight  segments  provided  the 
dominant  Mach  number  trends  are  contained 
in  the  wind  tunnel  data.  This  procedure 
may  require  significant  reprogramming 
from  vehicle  to  vehicle  depending  on  the 
particular  forms  of  the  wind  tunnel  models. 


Each  of  thes5  procedures,  as  well  as  certain  combinations 
of  them, should  be  investigated  using  appropriate  controlled 
test  cases.  Such  studies  are  recommended  for  future  effort. 


SECTION  VI 

AERODYNAMIC  MODEL  STRUCTURE  IDENTIFICATION 


BACKGROUND 

The  need  for  aerodynamic  structure  identification 
arises  when  the  form  of  the  missile  dynamic  model  is  uncertain. 
The  design  of  an  extended  Kalman  filter  to  estimate  aerodynamic 
coefficients  is  based  upon  a s’lecific  stochastic  model  for  the 
missile's  motion.  If  this  "filter  model"  is  inconsistent  with 
actual  missile  behavior,  poor  aerodynamic  coefficient  estima^ 
tion  accuracy  may  result.  The  objective  of  structure  identifi- 
cation is  to  determine  which  model  among  a given  class  of  models 
best  represents  the  physical  system  of  interest.  Thus,  the 
structure  identification  problem  is  one  of  approximation;  i.e., 
determining  which  model  within  the  given  class  best  approxi- 
mates the  input-output  response  of  the  physical  system.  A 
method  has  been  developed  for  identifying  model  structure  when 
its  uncertainty  consists  of  several  different  models  which 
might  be  best  or  correct.  This  method,  based  upon  hypothesis 
testing  theory,  is  demonstrated  and  evaluated  in  this  chapter. 
Further  examples  of  structure  identification  using  this  method 
are  given  in  (Reference  16). 

Structure  identification  is  an  important  component  of 
system  identification  and  parameter  identification  for  the 
following  reasons; 

• If  the  chosen  structure  is  too  complex 
(i.e.,  if  too  many  parameters  are  in- 
cluded in  the  model),  estimates  of  the 
relevant  parameters  (i.e.,  those  param- 
eters needed  to  described  the  system 
response)  will  be  degraded.  This  par- 


tjrulur  dil'lii'uHy  scvorc  when 

tlUTO  is  a limited  amount  oi  i ninil 
data  available 

• li  the  strueturo  ehosen  is  too  simple 
or  incorrect , the  physical  system  can- 
not be  adequately  represented  lor  any 
values  ol  the  parameters  in  the  selected 
mode  1 . 


Researchers  in  the  field  of  structure  identification 
have  recoitnixod  the  obvic'us  trade-('iff  that  exists  between  the 
number  of  parameters  used  in  a model,  and  the  ability  of  that 
mc>dol  tc'  fit  a ttiven  set  of  data  (. Rt'f erences  17-10).  Many  of 
the  techniqties  currently  in  use  for  selecting:  a model  are  based 
upon  least  sqtmres  fits  of  mcidel  parameters  to  the  data.  The 
model  structure  is  chosen  according  to  subjecTive  stuiistical 
criteria  which  depend  upon  variables,  such  as  "risk  levels," 
that  are  specified  in  an  ad  hoc  fashion.  In  addition,  least 
squares  approaches  are  theoretically  deficient  in  that  they 
are  based  upon  certain  assumptions  about  error  statistics  that 
are  rarely  satisfied  i n practical  problems  of  structure  iden- 
tification. When  these  assumptions  are  not  satisfied,  the 
least  squares  estimates  develop  bias  errors  which  can  severely 
hinder  the  model  selection  process. 

The  motivation  for  the  technique  used  here  comes  from 
previous  work  done  on  the  parameter  identification  part  of  the 
system  identification  problem.  Least  squares  techniques  were 
initially  employed  for  parameter  identification  due  to  their 
relative  simjil  icity ; however,  they  often  produce  biased,  in- 
accurate estimates.  Subsequently,  more  accurate  techniques 
evolved  such  as  maximvmi  likelihood  fe.g.,  as  applied  in  Ref- 
erence 20),  and  the  extended  Kalman  filter  (References  2 and 
12).  Inaccurate  results  have  also  been  observed  from  applica- 
tion of  least  squares  procedures  to  the  structure  identifica- 
tion problem  (References  21  and  22);  specifically  there  is  a 
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tendency  to  select  models  having  too  many  parameters.  In  this 
chapter  a new  technique  based  upon  extended  Kalman  filtering 
is  presented.  It  is  used  to  obtain  an  improved  solution  to 
the  problem  of  structure  identification  applied  to  missile 
aerodynamics . 

The  structure  identification  procedure  used  in  this 
study  is  illustrated  in  Figure  24.  A number  of  hypotneses , 

H. , , ...  H , are  defined,  where  H.  is  the  hypotheses  that 

± ^ X\  i , 

the  model  structure  upon  which  the  i''^  filter  is  based  is  best. 
The  choice  of  the  set  of  hypotheses  (or  set  of  candidate  mod- 
els) is  based  upon  the  a priori  information  available  about 
the  system  to  be  identified.  The  best  model  structure  can  be 
identified  by  operating  n extended  Kalman  filters  in  parallel, 
as  indicated  in  Figure  24,  with  each  filter  designed  according 
to  one  of  tne  candidate  models.  At  the  same  time  it  is  pos- 
sible  to  recursively  compute  the  probability  that  the  i 
model  is  correct;  making  use  of  the  filter  state  estimates 
X.  and  l.he  measurement  data;  the  relative  magnitudes  of  these 
probabilities  provide  a basi.s  for  selecting  the  model  which 
best  represents  the  data.  The  mathematical  details  of  this 
procedure  are  presented  in  Appendix  B. 

The  remainder  of  this  section  presents  a structure 
identification  study  for  an  air-to-air  missile.  Several  can- 
didate models  are  selected,  having  different  sets  of  aerody- 
na.mic  parameters.  One  of  these  models,  the  "truth  model",  is 
used  to  generate  synthetic  measurement  data.  Extended  Kalman 
filters  are  then  constructed  based  upon  each  of  the  candidate 
models  and  used  to  process  the  measurement  data.  The  prob- 
ability that  each  model  (each  hypothesis)  Is  the  correct  one 
is  computed  using  the  structure  identification  algorithm. 

The  ability  of  this  procedure  to  detei'mine  the  correct  system 
model  is  demonstrated  by  a number  of  simulations. 


M£ ARUREMCNI 
DATA 


Figure  24. 


Functional  Diagram  of  the  Structure 
Identification  Algorithm 


EXPERIMENTAL  RESULTS 


The  experimental  results  presented  in  this  section 
are  based  upon  the  air-to-air  missile  model  described  in  Sec- 
tion V.  Three  versions  of  this  model  are  employed  here  as 
truth  models  to  generate  synthetic  measurement  data.  The 
first  has  11  aerodynamic  coefficients,  the  second  has  16 
and  the  third  20.  Table  31  lists  the  paranieter  values  for 
the  three  truth  models.  The  control  inputs,  measurements, 
noise  sources,  dynamic  initial  conditions,  and  other  model 
properties  are  the  seme  as  the  model  fully  defined  in  Section 
V.  Note  that  Model  1 ha.-'  only  linear  aerodynamic  parameters, 
while  Model  2 has  five  additional  nonlinear  parameters,  and 
Model  3 has  four  more  nonlinear  parameters  than  Model  2. 
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TABLE  31.  TRUTH  MODEL  PARAMETER  VALUES 


I 

I 

I 


Before  proceeding  with  specific  test  cases  it  is  im- 
portant to  emphasize  the  difference  between  the  truth  model 
and  filter  mc'dels  in  a given  simulation.  The  truth  model  is 
used  to  generate  a specific  set  of  synthetic  measurements. 

The  filter  model  l.s  the  model  upon  which  a particular  filter 
design  is  based.  In  each  case,  for  a particular  measurement 
sequence,  two  or  more  extended  Kalman  filters  are  used  to  pro- 
cess the  data.  Each  of  these  filters  is  based  upon  one  of  the 
models  defined  in  Table  31.  The  outputs  of  these  filters  are 
then  compared  by  the  structure  identification  algorithm  to 
determine  which  filter  model  best  represents  the  truth  model 
used  to  generate  the  data.  The  initial  parameter  statistics 
for  the  filters  are  given  in  Table  32.  Other  than  these 
initial  statistics,  tne  filter  design  and  impleinentation  are 
the  same  as  the  nominal  case  described  in  Section  V. 

In  the  first  test  case  of  the  structure  identifica- 
tion procedure,  Model  1 is  used  to  generate  the  synthetic  mea- 
surement data.  Three  hypotheses,  and  are  formulated 

and  defined  with  their  corresponding  probabilities  as  follows: 

+ Vi 

H.  - Hypothesis  that  the  i filter  model  corre- 

^ spends  to  the  truth  model,  for  i = l,  2,  and  3 

P - Probability  that  the  hypothesis  H.  is  true, 
given  the  measurement  data.  ^ 

To  compute  the  three  probabilities  of  interest,  the  measure- 
ment data  are  processed  independently  by  three  extended  Kalman 
filters,  each  designed  in  accordance  with  a different  hypoth- 
esis (i.e.,  filter  model).  The  estimates  and  error  covariance 
matrix  of  each  filter  are  then  used  by  the  algorithm  described 
in  Appendix  B to  compute  the  three  hypothesis  probabilities  as 
a function  of  time. 
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TABLE  32.  FILTER  MODEL  PARAMETER  INITIAL  STATISTICS 


INITIAL 

UEAN 

INITIAL 

STANDARD  DEVIATION 

PARAUETER 

ALL  FILTER 
MODELS 

FILTER 
MODEL  1 

FILTER 
MODEL  2 

FILTER 
MODEL  3 

C 

’‘o 

-1.63 

0.4 

0.4 

0.4 

Cv. 

‘ a 

-•13.55 

10.0 

10.0 

10.0 

C,. 

* 5r 

-0.04-1 

0.01 

0.01 

C.Ol 

-4  3.. 55 

10.0 

10.0 

10.0 

-0.044 

0.03 

0.01 

0.01 

‘"M 

*^a 

•2G0. 3G 

50 . 0 

50.0 

50.  j 

2.9 

0.7 

0.7 

0.7 

-12500.0 

3000.0 

3000.0 

3000.0 

200.55 

50.0 

50.0 

50.0 

-2.9 

0.7 

0.7 

0.7 

C.N 

r 

-12500.0 

3000.0 

3000 . 0 

3000.0 

0.0 

- 

10.0 

10  0 

Cxg2 

0.0 

- 

10.0 

10.0 

3 

“ 

0.0 

~ 

1000.0 

1000.0 

0.0 

- 

.1300.0 

1000.0 

Cn,3 

0.0 

- 

1000.0 

lOOC.O 

0.0 

- 

- 

20n  . 0 

^11  ^ 
a Jr, 

0.0 

- 

I " 

^Na‘; 

0.0 

- 

- 

1000. 0 

^N-2. 

: vF 

0.0 

- 

- 

10.0 

Figure  25  illustrates  the  result  of  the  above  pro- 
cedure where  the  initial  a priori  probability  of  each  hypoth- 
esis is  assumed  to  be  one-third.  Note  that  initially,  the 
probabilities  remain  constant  and  equal  until  approximately 
one-half  second  into  the  flight.  During  this  period,  the 
missile  is  in  steady-state  flight  with  zero  control  inputs, 
(see  Figures  14  and  15).  This  behavior  indicates  that  the 
data  contains  little  or  no  information  that  can  be  used  to 
discriminate  between  the  three  hypotheses.  At  about  one-half 
second  the  control  input  begins  and  the  probability  of  the 

correct  hvpothosi.s,  P„  . rapidly  goes  to  1.0  as  the  others 

"l 

approach  zero. 

For  the  same  meusurerrient  data  as  above,  two  addi- 
tional cases  were  run  using  the  structure  identification  pro- 
cedure. Figure  26  illustrates  the  discrimination  achieved 
between  hypothesis  and  and  Figure  27  presents  a similar 
result  for  and  • In  each  case  the  correct  hypothesis  is 
rapidly  identified  soon  after  the  beginning  of  the  missile 
control  input. 

In  u second  case  study,  Model  2 is  used  to  generate 
the  measurement  data.  Figure  28  shows  the  results  of  the 
structure  identification  procedure  for  the  three  hypotheses. 
Note  that  the  algorithm  rapidly  identifies  the  correct  hypoth- 
esis, Ph3  remains  consistently  higher  than  because 

hypothesis  contains  all  of  the  parameters  of  Hg . whereas 
conta  .ns  only  a part  of  the  parameter  set  needed  for  Model  2. 

As  a further  illustration.  Figure  29  shows  the  re- 
sult of  computing  P^^  and  Pjig  from  the  Model  2 data.  Note, 
that  in  this  case  neither  hypothesis  is  correct.  Since, 
however,  filter  Model  3 contains  Model  2 as  a subset , whereas 
filter  Model  1 does  not,  the  probability  calculation  picks 
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Figure  29.  Probability  that  Hj  is  True  when 
H2  is  True,  for  i * 1 and  3 

That  hypothesis  which  best  represents  the  actual  data.  This 
demonstrates  that  the  technique  will  not  only  select  the  cor- 
rect hypothesis  if  it  exists,  but  it  tends  to  select  the  best 
or  most  likely  hypothesis  if  the  correct  hypothesis  is  not  known. 

In  the  final  test  case,  Model  3 is  used  to  generate 
the  measurements  and  the  sanie  three  hypotheses  are  tested. 

Figure  30  shows  that  again  the  correct  hypothesis  is  rapidly 
Identified. 

The  above  results  clearly  demonstrate  the  ability  of 
the  structure  identii icatlon  procedure  to  select  the  correct 
sy.stem  model  when  the  data  ccntain.s  informal ior.  that  allows 
the  discrimination  to  be  made.  In  practical  cases,  the  tech- 
nique can  be  used  to  evaluate  the  accuracy  oi  a series  of 
models,  none  of  which  may  .'luequnt cly  described  the  system 


Figure  30.  Probability  that  is  True  when 
Hj  is  True,  for  i = 1,  2,  and  3 

which  generates  the  data.  The  computat ioi^ial  cost  of  this 
procedure  is  simply  the  cost  of  applying  an  EKF  to  the  data 
for  each  proposed  model,  since  the  probability  calculations 
require  little  additional  computer  effort. 


SECTION  VII 

GBU-15  FLIGHT  TEST  DATA  PROCESSING 


BACKGROUND 

The  GBU-15  is  a family  of  nonthrusting  air-to-ground 
standoff  weapons  launched  from  an  aircraft  and  steered  during 
flight  so  as  to  impact  a designated  target.  Among  the  GBU-15 
family  of  modular  glide  weapons,  the  PWV  (Planar  Wing  Weapon) 
is  characterized  by  a high  aspect  ratio  swept  wing,  cruciform 
tail  and  bottom  mounted  fin.  Figure  31  gives  the  signifi- 
cant physical  characteristics  of  this  airframe.  Extensive 
wind  tunnel  testing  (over  250  hours)  has  been  conducted  in 
order  to  determine  the  aerodynamic  model  for  autopilot  de- 
sign and  performance  evaluation.  In  addition,  flight  tests 
have  been  conducted  to  exercise  the  airframe  throughout  the 
flight  regime  and  to  verify  performance  predictions. 

The  data  required  for  modeling  the  GBU-15  was  pro- 
vided by  AFATL  and  included  mass,  CG  and  inertia  characteris- 
tics, and  the  detailed  wind  tunnel  aerodynamic  model.  This 
wind  tunnel  model  includes  over  60  static  and  dynamic  terms, 
each  of  which  requires  a one-,  two- , or  three-dimensional  table 
lookup.  In  spite  of  the  complexity  of  this  model,  compari- 
sons of  in-flight  measurements  with  computer  simulated  mea- 
surements show  considerable  discrepancy  - particularly  in 
regions  where  the  vehicle  exhibits  rapid  dynamic  motions. 
Several  configuration  modifications  were  made  after  wind 
tunnel  testing  which  were  not  modeled  by  the  wind  tunnel  data 
a telemetry  antenna  was  located  aft  of  the  bottom  mounted  fin 
and  an  a-e  boom  was  mounted  in  the  nose  of  the  vehicle. 


AIRFOIL  SECTION  • NACA  65-4)0 
THEORETICAL  WING  AREA  - 16.59  FT^ 
WING  SP.AN  - 136.3  INCHES 
WING  CHORD  • 18.48  IN.  (STREAMWISEI 
WING  CHORD  • 16.00  (CHOROWISEI 
SWEEP  ANGLE  • 30  DEG 
ASPECT  RATIO  - 7.74 
INCIUENCE  ANGLE  • 2.6  (STREAMWISEI 
3.0  (CHOROWISEI 
OIHF.ORAL  ANGLE  • -10  OEG 
OVERALL  length  - 164.60  INCHES 
TAIL  SPAN  - 66  INCHES 
SWEEP  ANGLE  « 25  DEG 


136.30 


Figure  31.  Physical  Characteristics  of  the  GBU-15 


FLIGHT  TEST  DESCRIPTION 

All  data  from  a single  GBU-15  flight  test  was  pro- 
vided for  use  with  the  parameter  identification  algorithm 
discussed  previously.  This  flight  test  (AFATL  designation 
PPV-4)  consists  of  a high  altitude  launch,  followed  by  an 
altitude  hold  segment,  a maximum  L/D  glide  segment,  and  a 
rapid  pitch  down  to  a near  vertical  impact.  Several  roll 
and  yaw  maneuvers  were  executed  during  the  flight  to  verify 
autopilot  responses.  In  addition,  two  significant  pitch 
maneuvers  occurred,  one  at  the  transition  from  altitude  hold 
. to  glide  and  another  at  the  pitch  down  to  irapacc . 

i 


I 


The  I'liKht  test  data  consists  of, 

• A data  tape  containing  dipit izod  telemetry 
records  from  onboard  instrumentation  at 
approximately  a 20  hertz  data  rate 

• Atmospheric  data  from  instruments 
launched  prior  to  flight  test.  This 
data  includes  wind,  temperature  and 
atmospheric  density  measurements  as 
functions  of  altitude 

• Radar  data  processed  by  an  AFATL  software 
package. 

The  G0L:-15  onboard  instrumentation  included  a 
lateral  acceleronieter  and  normal  accelerometer,  three  body- 
fixed  rate  gyros,  a free  gyro  indicating  inertial  roll  and 
heading,  an  boom,  and  control  surface  deflection  sensors. 
Information  on  accelerometer  accuracy  and  a-S  boom  wind 
tunnel  calibration  was  provided. 

The  smoothed  radar  data  was  produced  by  an  AFATL 
missile  test  range  data  processing  package  which  utilized 
the  atmospheric  data,  together  with  various  radar  inputs, 
to  produce  estimates  of  missile  trajectory  characteristics. 
Fixed  interval  smoothing  methods  were  used  to  provide  esti- 
mates of  position  and  velocity  with  respect  to  inertial 
space  and  with  respect  to  the  relative  air  mass. 

The  telemetry  data  was  plotted  and  visually  in- 
spected in  order  to  determine  segments  suitable  fo:-  aero- 
dynamic parameter  identification.  Three  segmcui."'  , ruch  f 
about  ten-sec  length,  were  selected  based  on  l.<.  Jynair.ie 
activity  occurring  during  the  trajectory.  SeK  v.  ed  tele- 
metered flight  measurements  are  given  in  Appendix  C for 
each  of  the  three  segments.  Segment  I includes  a large 
pitch  maneuver  at  the  transition  to  constant  glide,  Seg- 
mient  11  includes  a yaw  maneuver  with  roll  coupling,  and 
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Segment  III  occurs  during  the  pitch  down  to  impact  and  in- 
cludes a large  pitch  plane  maneuver  with  some  roll  and  yaw 
activity.  Table  33  summarizes  the  initial  and  final  tra- 
jectory characteristics  from  the  radar  data  for  each  of  the 
three  segments. 

The  wind  magnitude  and  heading  from  the  atmospheric 
data  are  given  in  Figure  32.  The  three  flight  segments  occur 
at  altitudes  from  15,000  to  24,000  feet.  During  this  period 
the  wind  heading  and  magnitude  are  nearly  constant  with  alti- 
tude so  that  a constant  wind  magnitude  of  50  knots  with  head- 
ing of  305  degrees  east  of  north  was  assumed.  The  missile 
initial  heading  was  326.6  degrees  east  of  north  and  radar 
data  indicates  the  heading  was  maintained  closely  throughout 
the  flight.  Each  flight  segment  was  initialized  with  zero 
initial  heading,  and  with  inertial  wind  approaching  21,6 
• degrees  from  the  left  at  50  knots.  At  nominal  missile  veloci- 
ties, this  lateral  wind  component  can  contribute  about  3 
degrees  to  the  sideslip  (whereas  peak  measured  6 values  are 
about  5 degrees).  Thus,  wind  modeling  should  play  a key  roll 
in  the  yaw  plane  representation. 


TABLE  33.  SUMMARY  OF  RADAR  DATA  PROCESSING 
FOR  THREE  FLIGHT  SEGMENTS 


TRA.IFCTitin 

VARMBI.r 

StGMFNT  I 

.SFOMENT  J I 

SrCMENT  III  j 

v-io 

t-P 

t-30 

t“0 

t-Ht 

A) I i t udr  ( f t ) 

23f18g 

21  752 

19985 

1B14R 

19160 

15609 

Inrrt  lal  Vrlorlt v 

Maqn 1 1 udo  ( f ps  ) 

rj5o 

G2P 

072 

GR5 

4 78 

64  2 

R“1 nl j vr  Vcloc 1 1 V 

Maptiitudc  (fps) 

623 

002 

747 

767 

552 

676 

Mnch  Niimbor 

P.601 

0.659 

0 . 706 

0. 71!i 

0.020 

vi.nr.T 

Flight  Path 

Anplr  (dfp) 

-5.8-1 

-30.426 
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The  measured  atmospheric  temperature  was  used  to 
compute  the  actual  values  tor  speed  of  sound.  A comparison 
of  these  computed  values  and  the  standard  atmosphere  values 
are  given  in  Figure  33,  as  well  a.s  a comparison  of  measured 
and  standard  air  density.  The  standard  values  were  used  in 
the  identification  algorithm. 
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Figure  32. 


Measured  Wind  Data  for  GBU-15  Test  Flight 
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Figure  33.  Compurison  ol  Measured  Speed  of  Sound 
and  Density  With  Standard  Valutas 
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The  GBU-15  mas.'-  properties  consist  of  mass,  CG 
location,  and  inertias.  These  data  are  given  in  Table  34. 

The  lateral  and  normal  uccel erometei's  arc  located  near  the 

* 

rear  of  the  weapon  at  an  average  location  as  shown  in  the 


’♦The  two  aroeleromtiters  are  actually  located  about  9 inches 
apart  along  the  x-axis. 


TABLE  34.  GBU-15  AIRFRAME  DATA 


REQtIRED  AIRFRAME  DATA 

GBt-15  VALUE 

Mass  (slugs) 

86,354 

x,y,2  CG  Location  (in) 

0.875,  0.0,  -1.54 

O 

Roll  liiertia  islug-ft**) 

119.38 

Pitch  Inertia  (slug-ft'^) 

675.55 

2 

Vaw  Inertia  (slug-ft  ) 

747.42 

x,y,2  Accelerometer  Location (in) 

-55.93,  1.15,  4.65 

table.  Instrument  and  CG  location.s  are  given  with  respect 
to  the  aerodynamic  reference  point  in  body  axes,  defined  in 
Figure  31. 

The  a priori  aerodynamic  model  for  the  GBU  requires 
a parameterization  of  the  detailed  tabular  data  provided  by 
AFATL,  as  well  as  the  assignment  of  uncertainties  to  the  re- 
sulting parameter  values.  This  simplified  model  was  con- 
structed based  on  discussions  with  aerodynamicist s familiar 
with  the  GBU-15,  previous  e.xperience  with  winged  aerodynamic 
vehicles,  inspection  of  the  actual  flight  data  and  analysis 
of  the  wind  tunnel  data. 

A simplified  EKF  model  was  used  for  the  initial 
trials,  assuming  pitch  aerodynamics  to  Ve  decoupled  from  the 
yaw  and  roll  aerodynamics.  Yaw  and  roll  are  known  to  be 
highly  coupled  so  that  all  linear  yaw-roll  aerodynamic  ef- 
fects frequently  encountered  in  aircraft  are  modeled  (Ref- 
erence 23).  Additional  parameters  were  ’ -'vest  igated  based 
on  the  initial  identification  trials . 

Figure  3'  shows  representative  examples  of  the  wind 
tunnel  data.  Figures  34(a)  and  34(g)  demonstrate  the  complex 

i 
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(a)  Axial  Force  Coefficient 
vs  Angle-of-Attack 


(c)  Yaw  Moment  Coefficient 
vs  Sideslip 


Figure  34 . 
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(b)  Normal  Force  Coefficient 
vs  Angle-of-Attack 


(d)  Pitching  Moment  Coefficient 
vs  Angle-of-Attack 
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(g)  Roll  Moment  vs  Angle-of- 
Attack  and  Sideslip 


Figure  34. 


Representative  Wind  Tunnel  Data  (Concluded) 


nature  of  the  aerodynamics  for  this  vehicle.  Only  slight 
coefficient  variation  with  Mach  number  below  0.75  is  evident, 
so  that  coefficient  invariance  with  Mach  number  is  reasonable 
for  Mach  number  ranges  given  in  Table  34,  The  parameters 
investigated  in  this  study  are  given  in  Table  35  along  with 
the  initial  values  and  uncertainties. 


The  only  measurements  used  for  identification  during 
this  study  are  the  two  accelerometers,  three  rate  gyros  and 
the  a-5  sensors.  A conservative  initial  estimate  of  measure- 
ment noise  was  obtained  by  observation  of  the  actual  data. 
These  initial  estimates  were  modified  by  inspection  of  the 

EKF  residual  process.  The  accelerometer,  rate  gyro  and  a-B 

2 

sensor  noises  were  initially  given  by  2 ft/sec  , 0.01  rad/sec 
and  0.01  radian,  respectively. 


The  force  perturbation  inputs  (process  noise)  were 

2 

assigned  values  of  2 ft/sec  . The  moment  perturbation  inputs 
were  given  values  representative  of  a 0.3-degree  noise  in  con- 
trol deflection  measurement  at  a dynamic  pressure  of  250 
2 

Ib/fx  . Ail  studies  in  this  section  assumed  a wide-band 
(white)  process  noise  input. 


The  transfer  functions  relating  a and  6 to  the  con- 
trol inputs  6^  and  were  computed  from  the  mass  properties 
and  the  aerodynamic  coefficients.  This  ’’open  loop  airframe 
response"  indicates  lightly  damped  pitch  and  yaw  plane  re- 
sponses with  damping  ratios  of  approximately  0.18  and  0.08, 
respectively,  and  a natural  frequency  of  about  2.7  rad/sec  in 
both  planes.  Step  responses  from  the  6-DOF  a priori  airframe 
model  correlated  well  with  these  computations. 


The  GBU-15  utilizes  both  accelerometer  and  pitch 
rate  feedback  for  stability  agumentation.  Well  damped 


I 

TABLE  35.  GBU-15  PARAMETER  DATA  i 


•All  angles  are  expressed  in  radians. 


responses  to  pitch  and  yaw  plane  acceleration  conimand.s  are 
observed  with  .settling  times  of  about  1.5  and  2.5  seconds, 
respectively.  Some  concerns  have  been  expressed  in  the 
literature  about  the  ident if iability  of  dynamic  systems  im- 
bedded in  closed  loop  systems  (Reference  24),  although  studies 
discussed  in  Reference  24  indicate  successful  identification 
in  systems  operating  in  both  open  loop  and  closed  loop  modes. 


No  specific  investigations  into  the  effects  of  the  closed 
loops  on  identification  of  the  GBU-15  was  made;  all  measure- 
ments result  from  the  operational  closed  loop  system. 

IDENTIFICATION  RESULTS 

Initial  identification  tests  were  performed  using 
only  ten  measurement  samples  order  to  isolate  extreme  a 
priori  modeling  problems.  During  this  stage  the  data  rate 
of  10  hertz  (every  other  data  point)  was  selected,  dynamic 
state  initial  conditions  of  each  data  segment  were  adjusted 
and  measurement  noise  levels  were  modified  slightly  from  the 
a priori  values.  The  measurement  levels  assumed  in  the  EKF 
model  are  given  in  Table  36.  After  these  initial  trials,  a 
100-sample  identification  test  was  performed  for  each  of  the 
three  segments  discussed  previously.  The  21  parameters  in- 
dicated on  Table  35  were  identified  for  each  segment. 


TABLE  36.  ADJUSTED  MEASUREMENT  NOISE  LEVELS 


The  predicted  and  actual  measurement  comparisons, 
the  residual  processes,  and  the  predicted  residual  bounds 
which  occurred  during  parameter  identification  are  shown  in 
Figures  35  through  37.  The  predicted/actual  measurement  com- 
parison as  shown  in  these  figures  must  be  interpreted  care- 
fully. If  sufficient  process  noise  is  injected  into  a dy- 
namic model,  almost  any  set  of  parameter  values  may  produce 
an  excellent  predicted  measurement/actual  measurement  com- 
parison. In  addition,  the  dynaniics  which  generated  these 
measurements  utilized  time-varying  parameters  which  resulted 
during  the  course  of  parameter  identification.  A mo-^e  con- 
clusive demonstration  of  the  ability  of  the  final  parameter 
set  to  represent  the  data  is  given  later.  The  final  param- 
eter estimate  ana  the  predicted  standard  deviations  are  shown 
in  Table  37. 


Because  of  the  large  angle-of-attach  maneuvers  during 
data  segments  t and  III,  the  pitch  plane  aerodynamics  should 
be  highly  observable.  The  pitch  plane  dynamics  appear  to  be 
modeled  accurately  as  indicated  by  the  normal  acceleration, 
pitch  rate  and  angle-of -attack  residuals. 

Large  values  of  the  correlation  coefficients  between 
some  parameters  were  observed,  particularly  between  and 

Cz^  and  between  and  Segments  I and  III  p’roduoed 

correlation  coefficients  of  0.91  and  0.B7,  respectively, 
between  Cz^  and  High  correlations  such  as  these  make 

interpretation  of  actual  parameter  values  difficult  since  a 
family  of  parameter  values  may  produce  near  identical  residua^ 
processes.  The  estimate  of  was  consistent  between  the 

two  segments  and  showed  a slightly  positive  value  — indj eat- 
ing a less  stable  airframe  than  predicted  by  wind  tunnel  data. 
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Figure  35. 


Predicted  and  Actual  Measurement 
Comparison  for  Segment  I (Continued) 
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figure  3^.  Predicted  and  Actual  Measurement 

Comps risoi'i  for  Segment  II  (Continued) 
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TABLE  37 


FINAL  GBU-15  PARAMETER  E3T [MATES 
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Data  sogments  II  and  III  shoulO  provide  good  param- 
eter observability  for  the  coupled  roll  and  yaw  planes.  The 
residual  processes  from  these  segments  appear  to  indicate 
modeling  errors.  Note,  in  particular,  the  low  frequency  com- 
ponent exhibited  in  the  B residual  for  segment  II  in  Figure 
36(g).  A variety  of  attempts  were  made  to  improve  the  roll- 
yaw  model.  The  additional  coupling  terms  indicated  in  Table 
38  were  investigated  but  little  improvement  was  found.  Early 
in  the  study  some  confusion  arose  in  the  correct  signs  of  the 
input  control  deflections,  particularly  6r  and  5p.  Switching 
signs  on  these  inputs  produced  no  improvement.  Figure  37(g) 
suggests  the  sign  of  the  6 sensor  may  be  incorrect.  Changing 
this  sign  resulted  in  much  better  behavior  of  the  B residual 
as  shown  in  Figure  38. 


Several  factors  may  contribute  to  the  modeling  dif- 
ficulties in  the  roil-yav.  dynamics, 


• The  model  assumed  here  may  net 
accurately  capture  roll-yaw  effects. 

The  roll  wind  tunnel,  data  of  Figure 
34(g)  supports  this  premise. 

• The  initial  parameter  values  may  not  be 
accurate  enough  to  allow  the  filter  to  con- 
verge to  the  correct  values. 

The  yaw-roll  maneuvers  are  of  such  low 
amplitude  and  frequency  as  to  provide 
little  .information  in  the  duta  about 
the  yaw-roll  coefficients. 

• Signs  ol  the  measurements  and/or  control 
deflections  may  be  wrong. 


It  is  not  readily  apparent  from  the  data  presented 
thus  far  that  the  models  resulting  from  parameter  identifi- 
cation offer  an  improvement  over  the  a priori  models.  In 
order  to  demonstrate  that  the  collection  of  final  parameter 


135 


dnssois 


T/»BLE  :3f?.  ADDITIONAL  PARAMETERS  INVESTIGATED 


ADltiriOHAL 

PAy4M£,rt3S 

JKA.'tSTnAttO 

INjlTIAO 

VAata 

UACEhTAINTy  (o) 

0 

6 

i 

% 

0 5 

‘'h 

0 

0.2 

0 

0.05 

or 

S 

pa 

-O.rDS 

- 

c-;j 

0 

0.2 

“Sp 

V 

:i 

0 

5-0  (^PC) 

v *■ 
y 

0 

50  (lp6) 

?;  ■' 
i. 

10  1 

TWW«W\MWr*«<*WMU 

"'i'itti'i  compc'oeats . 

IKK  aiwii.v  wkixwm 

O.’OO  T-'’ 


'7r»  *1 


0 0'/5  ■ 


V 


•0.025  -i 


•0  OiO 


4 0 


< 0 (-  K'  >? 

TjMf.  SP.CrMfiWT  m (w) 


on. 

5 

l 

I 

•'«  'S 


¥ § 

p<?>  'S' 


i-| 

•c  w 

.-V— ■-■— -wp.;..--..!*.., 


^ « 6 TO  1} 

TlVk  If.v-J  ii.-.WENT  III  IKCI 


jUre  38,  CoiTnpar i'^on  of  TO'ealc'ced  ftvics  Actual  Sideslip 
lor  Se^.ment  HI  With  Input:  Measurernt-nt  Aa~ 
suiDcd  TO  P.er-'"e>st-T5t  Negative  Sidor,lip 


estimates  represents  an  Improved  model,  tv^o  predicted  mea- 
surement profiles  are  presented — one  generated  from  the 
initial  parameter  values  and  a second  from  the  final  param- 
eter estimates.  The  in provement  in  the  agreement  between 
the  actual  measurement  and  predictea  inoasurements  is  a mea- 
sure of  the  knowledge  gained  through  system  identification. 

As  demonstrated  in  Section  V the  initial  condition 
errors  and  process  noise  inputs  prevent  comparisons  resulting 
from  a straightforward  integration  of  the  dynamics  using  the 
initial  and  final  parameter  estimates.  Several  of  the  results 
noted  above  suggest  that  process  noise  plays  a significant 
role  in  the  models  which  have  been  fit  to  the  data.  A method 
of  comparison  which  treats  unknown  initial  conditions  and  prO“ 
cess  noise  consists  of  processing  the  measurement  data  via 
the  EKF  with  the  parameters  assumed  known,  i.e.,  no  parameter 
vector  is  augmented  to  the  dyn.amic  states.  If  the  EKF  is 
exc  cised  in  this  manner,  first  with  initial  parameter  values 
and  then  with  final  parameters,  two  measurement  profiles  are 
obtained  which  can  be  compared  with  the  actual  measurement 
data  to  show  system  model  improvement . 

The  complete  set  of  predicted  measurements  and  re- 
siduals from  data  Segment  I generated  by  the  EKF  using  param- 
eter estimates  fixed  at  their  initial  values,  is  shown  by 
Figure  39.  These  figures  indicate  considerable  residual  time- 
correlation  indicating  definite  system  modeling  errors.  The 
predicted  measurements  with  fixed  final  parameter  estimates 
are  given  by  Figure  40.  The  excellent  improvement  in  the 
residual  processes  demonstrates  the  increased  modeling  accu- 
racy achieved  with  the  parameter  estimates  obtained  by  the 
parameter  identification  algorithm.  The  initial  and  final 
parameter  values  from  Segments  II  and  III  were  not  compared 
in  this  manner,  though  it  is  expected  that  a siir.ilar  modeling 
improvement  would  be  indicated. 


137 


t*CRMAL  ACallEf^ATlOH 


YAW  RATE  |r«d/Mel 


(e)  Yaw  Rate 


(f)  Angle-of-Attack 


Figure  39.  Predicted  and  Actual  Measurement 
Comparison  With  Fixed  Initial 
Parameter  Values  (Continued) 


( ki ) Si  desl  i p 


Flpuro  liil.  PreciJcttni  :ind  Actual  Measurement 
Comp uri son  With  Fixed  Initial 
!’  a r amc  t e r V u 1 ut' s ( Co n c i u d c d ) 


EflATHX  UlTERAt. /kCCEL'RATION  Ct.' 


PiTCM  AATe  (r»d/«4cl 


(c)  Roll  Rate 


(d)  Pitch  Rate 


Figure  40.  Predicted  and  Actual  Measurement 
Comparison  With  Fixed  Final 
Parameter  Values  (Continued) 


I 

J 


^ 0 oor  *1  » 


% O.Olft  -1, 


I* 

I • I I 


7JMv  INTO  SROMENT  I Ucl 


0 9 4 e « 10  13 

TIME  INTO  SEGMENT  I (mc) 


(ej  Yaw  Rate 


'S 

^ 004 


5 0 0?  -1 


■\ 


'r#1tOivT<C'  [ 


•OVO 

0 7 f. 


f 10  1> 


TIWE  xNTO  StOMKNT  I 'ate) 


" -0.03 

3 


« 2 « e « «0  12 

TIME  INTO  segment  I Iwcl 


(f)  Angle-of-Attack 


? I ? 


Figure  40. 


Predicted  and  Actual  Measurement 
Comparison  With  Fixed  Final 
Parameter  Values  (Continued) 


144 


SECTION  VIII 


CONCLUSIONS  AND  RECOMJJENDATIONS 

SUMMARY  AND  CONCLUSIONS 

An  extended  Kalman  filter  algorithm  for  estimating 
aerodynamic  parameters  from  missile  flight  data  has  been  de- 
veloped and  evaluated  using  both  simulated  and  actual  flight 
test  data.  The  algorithm  includes  a general  purpose  6-DOF 
missile  airframe  model  suitable  for  representing  a variety  of 
missile  configurations.  Airframe  modeling  includes  the  effects 
of : 


• Time-varying  thrust  profile 

• Thrust  variation  with  altitude 

• Time- varying  mass  properties 

• Standard  atmospheric  model 

• Wind  versus  altitude  profile 

• Thrust  misalignment 

• CG  offsets 

• An  extensive  generalized  aerodynamic 
parameterizat ion . 

The  extended  Kalman  filter  requires  inputs  from 
various  flight  data  sensors.  The  measurement  options  inves- 
tigated in  this  study  are: 

• Body-fixed  rate  gyros 

• Body-fixed  accelerometers 

9 Angle-of-attack  and  sideslip  sensors 


Position  measurement 


• Attitude  measurements. 

Algorithm  verification  studies  and  filter  sensitiv- 
ity studies  have  been  performed  using  synthetic  data  from  a 
thrusting  short  range  interceptor  airframe  model.  Filtering 
performance  variations  resulting  from  various  measurement 
sets,  measurement  modeling  errors  and  measurement  noise  levels 
have  been  investigated.  Aerodynamic  representations  consist- 
ing of  a low  order  ’’linear”  aerodynamic  model  as  well  as  more 
extensive  nonlinear  models  including  up  to  36  aerodynamic 
parameters  were  also  studied.  Conclusions  resulting  from  the 
verification  and  sensitivity  studies  are: 

• The  extended  Kalman  filter  is  w'^ll 
suited  as  an  aerodynamic  parameter 
identification  tool  for  large  scale 
airframe  modeling. 

• The  computational  burden  — both  in  core 
size  and  computation  time  — are  not  ex- 
cessive for  the  realistically  modeled 
example  discussed  here. 

An  additional  study  was  to  be  performed  using  synthetic 
flight  data  provided  by  the  Air  Force  Armament  Laboratory. 

This  data  was  produced  by  a six-degrees-of-f reedom  interceptor 
simulation  developed  independently  from  the  model  discussed 
here.  No  parameter  identification  was  performed  as  a result 
of  the  large  aerodynamic  coefficient  variations  with  Mach 
number  during  the  flight.  All  aerodynamic  coefficients  varied 
on  the  order  of  30  percent  per  second  both  during  boost  and 
during  coast  deceleration.  Several  recommended  procedures  for 
aerodynamic  parameter  identification  with  large  Mach  number 
variations  are  given  in  Section  V. 


A structure  tdentificatlon  algorithm  also  has  been  de- 
veloped and  evaluated  using  synthetic  measurenient  data.  This 
algorithm,  used  in  conjunction  with  the  parameter  identifica- 
tion algorithm,  can  select  that  model  from  a candidate  family 
of  aerodynamic  models  which  most  likely  produced  the  input 
measurements.  The  capability  of  the  structuie  identification 
procedure  is  demonstrated  by  test  cases  in  which  the  algorithm 
consistently  and  rapidly  selects  the  correct  model  structure 
from  three  aerodynamic  representations  --  an  ll~para.meter 
model,  a 16-parameter  model  , and  a 20-parameter  model. 


The  identification  algorithm  has  been  used  to  process 
actual  flight  data  from  an  aerodynaniically  controlled  glide 
weapon.  This  experiment  provided  an  operational  test  of  the 
software  and  exercised  all  phases  of  the  aerodynamic  modeling 
procedure.  The  following  tasks  were  performed  in  preparation 
for  parameter  identification: 


• Interfacing  with  the  flight  test  agency 
so  that  missile  modeling  data,  flight 
measurement  tapes,  radar  data,  etc.,  could 
be  obtained,  interpreted,  and  placed  in  the 
proper  format  required  by  the  identification 
software . 

• An  a priori  airframe  model  for  use  in 
parameter  identification  v,^as  developed. 
Considerable  simplification  of  the  com- 
plex tabular  aerodynamics  was  required. 

<»  The  flight  data  was  visually  inspected  to 

isolate  segments  suitable  for  identification. 
Measurement  noises,  process  noises  and  state 
initial  conditions  were  approximated  for 
three  flight  segments.  Flexible  data  editing 
and  plotting  capabilities  were  developed  to 
assist  in  this  process. 


Parameter  identification  of  a 21-pai'ameter  model  was 
performed  for  the  three  data  segments.  Tlie  rc.'='.ulting  residual 


process  from  the  pitch  plane  measurements  appeared  well  be- 
haved indicating  pitch  plane  modeling  to  be  adequate.  The 
roll  and  yaw  plane  measurement  residuals  indicated  model  de- 
ficiencies in  the  coupled  roll-yaw  dynamics.  Various  attempts 
to  add  additioiiai  parameters  failed  to  improve  the  residual 
behavior.  Possible  causes  of  these  identification  difficul- 
ties are: 


• Tne  a priori  model  derived  from  the 

wind  tunnel  data  is  too  simplistic.  Com- 
plex aerodynamic  behavior  can  be  seen  in 
the  wind  tunnel  aerodynamic  model,  -which 
was  neglecxed  in  the  a priori  model  used 
for  identification. 

« The  airframe  roll-yaw  maneuvers  during  the 
three  selected  flight  segments  may  be  too 
low  in  amplitude  and  frequency  to  excite 
the  airframe  dynamic  modes  sufficiently  to 
al  lew  separation  of  the  aerodynamic  ef  fects 
iroir. process  and  measurement  noises. 

e The  parameters  may  be  initialized  incorrectly. 
Large  initial  parameter  errors  together  with 
statistically  inconsistent  initial  uncertainty 
levels  may  not  allow  parameter  estimates  to 
converge  to  the  "correct”  values. 

• Sign  errors  in  control  deflection  inputs, 
and/or  measurements  may  be  present  on  the 
input  data. 


Irrespective  of  the  di'f f iculties  indicated  above, 
the  overall  GBU-15  modeling  exercise  was  encouraging.  A dem- 
onstration of  modeling  improvement  after  identification  was 
given  for  which  the  extended  Kalman  filter  with  parameters  fixed 
at  their  estimated  values  was  used  to  process  flight  n’easure- 
ments.  These  results  were  compared  to  a similar  case  where 
the  parameters  were  fixed  at  their  initial  values.  A distinct 
inodel  improvement  was  demenstrated  by  comparing  the  residual 
process  behavior  from  these  two  cases. 
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RECOMMENI‘ATIONS 

All  stuci.jes  to  date  indicate  that  the  EKF  parameter 
identification  algorithm,  together  with  the  general  purpose 
missile  airframe  model  and  the  structure  identification  soft- 
v;are,  is  a useful  tool  for  aerodyaamic  modeling  from  flight 
data.  Several  recommend at  ions  for  improvements  and  modifica- 
tion to  the  existing  software  are  given  below, 

• Develop  a well-documented,  operational 
software  package. 

• Modify  the  software  to  allow  identifica- 
tion ol  missile  airframes  during  flight 
segments  containing  large  Mach  number 
variations.  Specific  recormnendations 
are  given  in  Section  V. 

• Include  the  capability  to  estimate  addi- 
tional airframe  parameters  such  as  CG 
offsets,  thrust  misalignments,  sensor 
misalignments  etc. 

• Add  a rani,  om  wind  model  and  perfcrm  addi- 
tional studies  into  the  offecte  of  wind 
gusts . 

Further  study  into  01>U-15  parameter  identii  icaticn 
is  also  recommended,  A major  omission  in  the  GEU'-l.'i  data 
processing  effort  was  the  failure  to  p-i^tjperly  develop  and 
verify  the  a priori  model  used  for  parameter  ident li icat i on . 

A more  extensive  a priori  model  is  required  and  the  resulting 
model  should  be  verified  against  existing  AFATL  GBU-15  simu- 
lations. This  exercise  woiild  ensure  a thoro(5p;h  un<lerstanding 
of  the  significant  features  of  the  wind  tunnel  model  and  pro- 
vide confidence  in  the  overall  accurac!/  of  tne  a pritiri  model- 
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APPENDIX  A 


IDENTIFICATION  ALGORITHM  OUTPUT 


This  appendix  presents  selected  plots  of  significant 
variables  computed  during  the  parameter  identification  studies 
of  Section  V.  Five  groups  of  figures  are  included.  These 
figures  are: 


A-1  State  estimation  errors  and  predicted 
one  sigma  bounds  for  the  11-parameter 
nominal  test  case 

A-2  Parameter  estimation  errors  and  predicted 
one  sigma  bounds  for  the  11-parameter 
nominal  test  case 

A-3  Re.Jidual  process  and  residual  one  sigma 

bounds  from  the  11-parameter  nominal  test 
case 

A-4  Parameter  estimation  errors  and  predicted 

one  sigmia  bounds  from  the  36-parameter  test 
case 

A-5  Residual  process  and  residual  one  sigma 
bound  for  the  36-parani  ' test  case 


The  solid  lines  shown  on  each  figure  represent  tne 
one  sigma  bounds  as  predicted  by  the  EKF.  The  vertical  dashes 
give  the  state  error,  parameter  error  or  residual  which  oc- 
curred at  the  particular  time  during  filtering. 
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(a)  STATE  ERROR 


(b)  Vy  state  error 


(..)  Vj  STATE  ERROR 


(d)  state  error 


Figure  A-l  . 


State  Estimation  Errors  for  the 
Nominal  Test  Case 
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Figure  A-1 


State  Estimation  Errors  for  the 
Noniinal  Test  Case  (Concluded) 
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Parameter  Estimation  Errors  for 
the  Nominal  Test  Case  (Concluded) 
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Figui’e  A-3. 


Re.sidual  Processes  from  the 
Nominal  Test  Case  (Coneluaed) 
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Parameter  Errors  from  the 
36-Parameter  Case  (Continued) 
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Parameter  Errors  from  the 
36  Parameter  Case  (Continued) 
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APPENDIX  B 


STRUCTURE  IDENTIFICATION  USING  THE 
EXTENDED  KALMAN  FILTER 


The  structure  identification  technique  used  in  this 
study  is  based  upon  the  extended  Kalman  filter  and  hypothesis 
testing  theory.  A number  of  hypotheses,  H.,  , H,,,  ...  H , are 
defined;  each  are  associated  with  a candidate  model  for  the 
system  structure  to  be  identified.  H.  is  the  hypothesis  that 
the  i model  best  represents  the  system.  The  choice  of  the 
set  of  hypotheses  (or  set  of  candidate  model  structures)  is 
based  upon  the  a priori  information  available  about  the  system 
The  best  model  structure  is  th  t structure  that  best  appro.^:!- 
mates  the  input-output  response  of  the  system.  This  structure 
is  found  by  processing  the  measurement  data  independently  with 
n extended  Kalman  filters,  as  indicated  in  Figure  B-1,  with 
each  filter  designed  according  to  one  of  the  candidate  models 
or  hypotheses.  It  is  then  possible  to  recursively  compute  the 
probability  Pj|^  that  the  i^^  model  is  best  or  correct,  making 
use  of  the  filter  state  estimates,  their  filter  computed  co- 
variances,  and  the  measurement  data.  The  relative  magnitudes 
of  these  probabilities  provide  the  basis  for  selecting  the 
model  which  best  represents  the  data  in  a statistical  sense. 

The  notation  used  in  Figure  B-1  is  defined  as  follows 

't'  V) 

H.  Hypothesis  that  the  i'"  system  model 

^ is  correct 

P„  (k)  Probability  that  the  i^^^  system  model 
is  correct  given  the  measurements  up 
to  and  including 


araifttiHBi^iVA*' 


Figure  B-1 


Extended  Kalman  Filter  Approach 
to  Structure  Identification 
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The  system  output  measurement  at  time,  tj^ 

The  predicted  value  of  z,  based  upon  i he 
data  up  to,  but  not  including  z,  . for  the 
EKF  based  upon  the  i^  model 

The  one-step  predicted  measurement  residual 
at  the  step  for  the  EKF  associated  with 
the  i't^  model. 


^ Vs 

The  output  of  the  i^  extended  Kalman  filter  immediately 

state  estimate  based  upon 


prior  to  the  measurement  is  a 


the  first  k-1  measurements  and  knowledge  of  the  past  and  current 
system  input.  This  estimate  is  used  to  predict  the  value  of  the 

the  1 


measurement  using  the  nonlinear  system  rutput  equation  for 
. th 


model,  i e. 


The  one  step  prediction  residuals  for  the  i 
by: 


,th 
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EKF  are  then  given 
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and  the  residual  covariance  by: 

. .T 
— k ~k 


,i 


'k 


Cov 


. . .T  . 

P?  (“)  + kI 

k k k K 


where 


3h 

. 

k 3x 


■■■  [2.] 


-k 


^k  ik^'^ 


(B-3) 


(B-4) 


i i “t  h 

and  is  the  error  covariance  of  x^(-)  computed  by  the  i^ 

filter . 


The  objective  is  to  use  the  above  filter  outputs  and 
equations  to  compute  the  probability  of  each  hypothesis.  For 
notational  purposes > let; 
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Appling  Bayes’  rule,  we  can  write: 
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Hated  upon  thi  assumpt  iout;  that  the  system  which  generates  the 
dai a is  d-iven  only  by  known  inputs  and  white  noise  and  the 
meas-irenients  are  corrupted  only  by  2ero  mean  white  sequences, 
then  the  residuals  should  be  zero  mean  gaussian  sequences  if 
the  .'Kl  linearization  assumption  is  valid.  In  this  case, 
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M 

(2Ti)^idet  S^i- 


(B-7) 


where  M is  the  dimension  of  the  residual  vector  and  [det 

denotes  the  absolute  value  of  the  determinate  of  the  matrix  sf. 

k 

Equation  B-6  is  used  recursively  to  compute  the  pro- 
ability of  each  hypothesis  for  i = 1,  2 •••  n.  It  must  be 

A 

initialized  at  assumed  a priori  values  for  the  probability  of 
each  hypothesis,  i.e., 


jj  (0)  - P(H^|(no  measurement))  (B-8) 


for  i = 1,  2,  •••  ii.  Based  on  k measurements,  the  decision  pro- 
cedure involves  selecting  the  model  or  hypothesis  which  has 


the  largest  probability. 


Ph.  (k). 


Several  comments  can  be  made  about  the  above  prc.<cedure: 


• The  technique  evaluates  each  model  on  the 
ability  of  the  associated  EKF  xo  predict 
future  observations  (i.e.,  it  is  based  up- 
on the  one-step  prediction  residuals) 

• The  decision  procedure  is  totally  objective 


I 


ii 

i 

-I 


The  tii'Pthod  provides  a direct  measure  of  the 
confidence  levels  {Py  (k)}.  Furthermore, 

“i 

the  probability  calculations  are  recursive, 
so  that  considerable  computational  savings 
can  be  incurred  in  cases  where  the  algorithm 
converges  rapidly 

The  technique  treats  the  multiple  hypotheses 
problem  directly 


In  dealing  with  a large-scale  system,  one 
can  perform  preliminary  calculations  to 
obtain  insight  as  to  the  oraer  of  the  sys- 
tem and  thereby  reduce  the  number  of 
models  that  need  be  considered.  Even  wjth 
a large  number  of  models,  man}/  of  the 
P,.  (k)'s  will  quickly  converge  to  zero  so 

that  the  associated  models  can  be  eliminated 
from  consideration  imnjediately 


The  developmen;i:.  of  the  extended  Kalman  fil- 
ter for  the  i^*'^  mc:del  is  based  on  the  as- 
sumption that  one  can  effectively  lii'.carize 
the  nonlineav  functions  about  the  state 
estimates  {xj;}.  little  is  known  about  the 
range  of  validity  of  this  assumption,  liow- 
ever,  because  of  the  genci'ally  good  re.sults 
that  have  been  obtained  with  extendoa  Kalman 
filters,  vis-a-vis  other  methods  in  param- 


they  will  offer  a significant  improvement 
for  structure  identit icati on . 


APPENDIX  C 

TELEMETRY  DATA  FROM  A GBU-l‘j  FLIGHT  TELT 


The  figures  contained  in  thif-r  appendix  represent 
typical  teleiActry  records  used  as  input  measurciT'ents  for  the 
system  identification  exercise  performed  with  the  GHL-lo 
standoff  glide  weapon.  Data  from  Segments  1.  II  and  III  ns 
defined  in  Section  VII  are  included.  Selected  records  of 
y-  and  ^.-axis  accelerometers , pitch,  y.aw  and  roll  rate  gyros, 
the  a-B  .sensor.s,  and  composite  pitch,  yaw  and  roll  tail  do- 
flection-'  are  shown.  All  digitized  records  are  plotted  at 
19.438  hertz  as  provi.ded  by  AFATL  with  data  points  corrected 
by  s t a igh t 1 ines  , 
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Figure  C-l(a).  Segment  I y-Axis  Acceleration 


Figure  C-l(b). 


Segment  I r.-Axis  Acceleration 
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Figure  C-3(a).  Segment  III  y-Axis  Acceleration 


Figure  C-3(b)  . Segment  III  z-Axis  Acceleration 
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Segment  III  Yaw  Control  Deflection 
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